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The  Hexapod  Coordinate  Measuring  Machine  (HCMM)  is  composed  of  six 
telescoping  struts  with  laser  interferometric  displacement  measurement  and  a  single 
fixed-length  center  rod.  The  laser  interferometers  measure  the  length  changes  of  each 
strut.  However,  in  order  to  calculate  the  coordinate  of  the  probe,  it  is  necessary  to  know 
the  absolute  length  of  each  strut  and  the  distance  between  the  three  joints.  Thus  there  are 
nine  kinematic  parameters  to  be  determined  in  the  calibration  process.  The  nine 
parameters  are  the  three  base  lengths  and  the  initial  lengths  of  the  six  struts. 

This  research  presents  a  method  to  find  the  initial  strut  lengths  and  the  base  joint 
locations  without  adopting  any  external  artifacts  or  measurements.  During  calibration  the 
measuring  probe  is  moved  to  a  number  of  different  positions  within  the  measuring 
volume.  For  each  position,  the  center  rod  length  is  assumed  to  remain  constant.  The 
length  is  calculated  based  on  an  initial  guess  at  the  machine  kinematic  parameters.  A 
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non-linear  least  squares  problem  is  then  formulated  to  search  for  the  values  of  the  base 
locations  and  initial  strut  lengths  which  minimize  the  variation  in  the  computed  center  rod 
length. 

A  computer  simulation  of  the  self  calibration  process  has  been  implemented. 
Noise  content  in  each  laser  interferometric  measurement  signal  was  monitored  to 
determine  its  effect  on  the  quality  of  the  calibration  results.  It  was  found  that  the 
magnitude  of  calibration  uncertainty  is  proportional  to  the  magnitude  of  the  noise  level. 
The  simulation  tool  also  found  that  the  data  collected  from  some  specific  motions  is 
insufficient  for  calibration  use.  Thus,  the  observation  indices  were  used  to  determine  the 
usability  of  a  data  set.  This  research  also  investigated  the  influence  of  the  number  of 
position  used  for  calibration  and  how  it  related  to  the  uncertainty  of  the  identified 
parameters.  Experiments  were  performed  to  estimate  the  thermal  time  constants  of  the 
HCMM  components.  Thus,  one  can  determine  the  time  period  that  the  HCMM  has  to  be 
re-calibrated.  In  the  final  of  this  research,  we  demonstrated  how  to  evaluate  the  HCMM 
volumetric  error. 


CHAPTER  1 
INTRODUCTION 

A  conceptual  five  degree  of  freedom  hexapod  coordinate  measuring  machine 
(HCMM)  was  proposed  by  Jokiel  and  Ziegert  [1].  The  research  presented  in  this  paper  is 
directed  toward  to  the  development  of  the  self-calibration  algorithm  for  the  HCMM. 

1.1  Conventional  Coordinate  Measuring  Machine 

Based  on  a  machine's  accuracy,  price,  work  volume,  and  throughput,  Slocum  [2] 
categorized  the  coordinate  measure  machines  (CMMs)  into  three  market  brackets: 
ultraprecision,  throughput,  and  small  manual  machines.  Those  CMMs  in  the 
ultraprecision  group  have  to  achieve  submicron  accuracy.  Only  a  few  machines  were 
being  sold  worldwide  because  of  their  high  cost  and  the  great  amount  hand  finishing 
involved  in  their  manufacture.  The  throughput  group  was  responsible  for  the  sale  of  a 
great  number  of  machines  worldwide.  These  machines  typically  display  accuracy  on  the 
order  of  400  to  600u.  inch  (10  to  15|im)  and  were  capable  of  making  rapid  computer 
controlled  measurements.  The  last  bracket  of  CMMs,  in  general,  was  made  up  of  small 
and  benchtop  type  machines.  The  market  for  these  smaller  machines  is  fairly  limited  and 
highly  competitive. 

Typically  a  CMM  uses  serial  kinematic  linkages,  each  providing  a  degree  of 
freedom.  The  configuration  could  be  a  bridge  type  or  cantilever  construction.  Generally 
these  designs  employ  a  three-dimensional  rectilinear  Cartesian  coordinate  system  since 
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the  design  is  based  on  three  mutually  perpendicular  prismatic  joints.  A  probe  system 
which  establishes  the  relationship  between  the  workpiece  and  the  machine  coordinate 
system  is  attached  at  the  end  of  last  serial  link.  The  commonly  used  bridge  type  CMM  is 
shown  in  Figure  1  - 1 . 


When  a  probe  touches  a  workpiece,  its  position  is  actually  determined  from 
displacement  transducers  which  are  embedded  along  each  of  the  three  axes  of  CMM. 
Thus,  the  structure  must  be  stiff  enough  to  avoid  deflections.  The  required  stiffness  is 
commonly  achieved  using  a  large  structure  with  high  mass.  These  large  masses  and  their 
associated  high  moment  of  inertia  must  be  driven  by  the  machine  drive  and  control.  A 
high-performance  CMM  is  characterized  not  only  by  accuracy,  but  also  by  the  ability  to 


Figure  1-1  Bridge  type  CMM 
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move  quickly  from  point  to  point  while  making  "measurements  on  the  fly."  Heat 
generation  from  such  a  high  speed,  high  power  driving  system  required  to  quickly  move 
large,  high  stiffness  machine  elements,  would  cause  significant  distortions  of  the 
machine.  In  order  to  provide  high  accuracy  measurements  the  CMM  has  to  be  placed  in  a 
temperature-controlled  environment  to  minimize  the  structural  distortion  due  to  thermal 
gradient  in  the  machine  elements.  Therefore,  the  total  cost  for  the  machine  itself  and  the 
creation  and  maintenance  of  the  required  working  environment  is  a  significant  economic 
issue  for  CMM  users. 

1 .2  Hexapod  Coordinate  Measuring  Machine 

An  alternative  is  the  hexapod  coordinate  measuring  machine  (HCMM).  This 
machine  uses  a  parallel  kinematic  structure  in  which  light  weight  telescoping  struts  are 
used  to  reduce  the  heavy  driving  power  needed  in  the  conventional  CMMs.  The  HCMM 
is  designed  to  self-calibrate  and  compensate  for  thermally  induced  changes  in  its  own 
elements.  Thus,  the  requirements  for  temperature  control  in  the  operating  environment 
are  greatly  reduced. 

The  HCMM  is  composed  of  six  telescoping  struts  and  a  single  fixed-length  center 
rod.  The  six  telescoping  struts  are  mounted  in  pairs,  with  one  pair  at  each  of  the  three 
spherical  joints  around  the  base  of  the  machine.  One  strut  from  each  pair  mounts  to  the 
top  of  the  center  rod,  and  one  strut  from  each  pair  mounts  to  the  bottom  of  the  center  rod. 
Each  telescoping  strut  carries  a  laser  interferometric  sensor  to  measure  its  length.  The 
configuration  of  the  machine  resembles  two  tetrahedrons  which  share  a  common  base  and 
whose  apex  points  are  connected  by  the  fixed  rod.  A  probe  is  attached  to  the  lower 


4 

sphere  such  that  the  probe  tip  is  collinear  with  the  center  rod.  This  kinematic  design 
gives  the  center  rod  five  degree  of  freedom  (5-DOF)  to  move  within  the  work  space.  The 
prototype  HCMM  is  shown  in  Figure  1-2.  The  expected  HCMM  system  performance 
was  outlined  by  Jokiel  in  his  thesis  [3]  as  follows 

1 .  Overall  accuracy:  ±89(1  to  ±213|Li  inch  (2.2  to  5.3  |im) 

2.  Workspace:  A  cylindrical  volume  15  inches  in  diameter  and  15  inches  in 
height. 

3.  Speed:  Maximum  constant  probe  speed  of  30  in/sec.  Strut  extension  speed  of 
25  in/sec. 

The  HCMM  can  be  divided  into  four  subsystems:  driving  and  control  subsystem, 
displacement  measurement  subsystem,  kinematic  (geometry)  subsystem,  and  calibration 
subsystem.  Each  of  the  subsystems  is  dedicated  to  achieving  a  critical  part  of  the  overall 
accuracy  and  performance  requirements.  A  brief  description  of  each  subsystems  follows: 
Driving  System 

A  prototype  servo  actuated  strut  of  the  HCMM  was  designed  and  manufactured  by 
Mevoli  [4].  The  maximum  and  minimum  strut  lengths  are  30  inches  minimum  and  50 
inches  maximum.  A  rotating  nut  type  drive  was  used  to  drive  a  hollow  ball  screw  and  to 
achieve  the  required  speed.  The  hollow  ball  screw  allows  the  laser  light  to  pass  through. 
The  strut  and  the  drive  system  are  shown  in  Figure  1-3.  Its  maximum  acceleration 
capacity  is  64.4  ft/s2  or  2G. 
Displacement  Measurement  System 

The  HCMM's  linear  displacement  measurement  system  is  modified  from  an 
instrument  called  the  Laser  Ball  Bar,  developed  by  Ziegert  and  Mize  [5].  It  consists  of  a 
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Figure  1  -2  The  hexapod  coordinate  measuring  machine 
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Figure  1-3  The  strut  and  actuating  system 


heterodyne  laser  interferometer  with  a  length  measurement  uncertainty  of  ±10|i  inch  or 
±0.25u\  meter.  Two  red  laser  waves  generated  by  HeNe  laser  with  nearly  the  same 
frequency  are  shot  out  from  an  emitter  which  is  located  at  one  end  of  the  strut  near  the 
base  sphere  to  a  moving  retroreflector.  The  moving  retroreflector  which  functions  as  a 
mirror  is  mounted  at  the  other  end  of  the  ball  screw  near  the  moving  sphere.  A  phase 
detector  is  used  to  measure  the  phase  difference  (A()))  between  the  laser  beam  as  it  leaves 
the  laser  head  and  the  reflected  beam.  Thus,  the  relative  distance  traveled  by  the 
retroreflector  is  found  as 


Ax  = 


A(p  X 
An 


(1-D 


where  A,  is  the  red  light  wave  length  (632.8  nm)  from  the  HeNe  laser.  Thus,  the  distance 
change  between  the  fixed  sphere  (or  base  sphere)  and  the  moving  sphere  could  be 
measured  by  the  above  system.  Figure  1-4  illustrates  how  the  heterodyne  laser 
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interferometer  works.  The  Af  1  shown  in  Figure  1-4  is  the  frequency  shift  caused  by 
Doppler  effect  when  the  retroreflector  moves.  The  advantage  of  the  laser  interferometer 
is  that  it  has  a  relatively  small  measurement  uncertainty  and  is  less  sensitive  to 
temperature  change  compared  with  other  displacement  measurement  devices.  However, 
it  has  to  be  noted  that  the  laser  interferometer  can  only  measure  relative  changes  in  the 
length  of  strut.  The  initial  strut  length  is  unknown.  Thus,  the  HCMM  needs  to  perform 
an  initialization  (or  calibration)  procedure  to  determine  the  initial  strut  length  before  it 
can  be  used  to  take  functional  measurements. 


fixed  cube  corner 
////////rj/. 
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Figure  1-4  Illustration  of  a  heterodyne  laser  interferometer 


Kinematic  System 

In  order  to  overcome  the  static  and  dynamic  deflections  under  actuation  forces, 
the  structure  has  to  be  very  stiff.  A  maximum  of  200  pound  force  acting  on  each  strut, 
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each  spherical  joint  and  each  of  the  three  supporting  columns  is  estimated  and  is  used  for 
determining  the  dimensions  and  materials  of  the  structure.  Angular  contact  bearings  are 
installed  on  all  the  spherical  joints  to  prevent  deflection  in  their  axial  direction. 

Because  the  spheres  in  the  base  and  moving  joints  provide  constant  references  for 
length  measurement,  the  imperfections  in  the  form  of  the  spheres  will  result  directly  in 
measurement  error  along  the  length  of  each  strut.  Therefore,  using  spheres  with  highly 
accurate  spherical  form  is  very  important.  For  instance,  a  grade  5  sphere  has  a  spherical 
form  tolerance  of  5.0fiin  (0. 125|im).  The  solid  triangle  frame  is  made  of  steel  to  provide 
a  heavy  thermal  mass  and,  consequently,  to  create  a  large  thermal  time  constant.  The 
symmetric  kinematic  structure  design  also  provides  symmetric  thermal  expansion  and 
thus  reduces  the  thermally  induced  distortion  of  the  structure. 

In  addition  to  the  measurement  accuracy  considerations,  many  other  factors  have 
to  be  considered  during  the  design  stages  of  the  HCMM  such  as  the  operation  safety,  part 
loading  convenience,  manufacturing  costs  and  the  ease  of  machine  installation. 
Calibration  System 

The  accuracy  of  the  HCMM  depends  on  the  knowledge  of  the  absolute  length  of 
each  side  of  the  dual  tetrahedron  kinematic  linkage.  If  there  is  any  measurement  error  on 
the  displacement  sensor  (or  the  laser  interferometer),  it  will  result  in  an  error  of  the 
measured  strut  length.  Therefore,  a  frequent  calibration  for  the  laser  interferometer  is 
important.  In  addition  to  the  first  level  calibration,  the  nominal  kinematic  parameters 
need  to  be  calibrated  frequently.  The  nominal  kinematic  values  may  be  obtained  from  the 
blue  print  of  the  design  drawing,  but  the  errors  from  manufacturing,  installation,  and 
thermal  effects  are  inevitable.  A  calibration  technique  for  finding  the  actual  kinematic 
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parameters  is  needed  before  the  machine  can  make  any  measurement.  The  development 
of  such  a  calibration  technique  is  the  primary  objective  of  this  research. 

1.3  Scope  of  Research 

As  mentioned  earlier,  the  HCMM  accuracy  depends  on  the  knowledge  of  the 
kinematic  parameters.  Though  the  laser  interferometer  measurement  system  is  used  in 
this  machine,  it  only  provides  the  relative  length  changes  of  the  struts.  The  initial  strut 
lengths  and  the  distances  between  base  spheres  must  be  determined. 

Finding  the  actual  kinematic  parameters  for  a  hexapod  machine  is  much  more 
difficult  than  for  a  traditional  Cartesian  type  machine.  Over  the  past  two-hundred  years 
people  have  established  a  mature  technique  to  find  the  parametric  errors  of  the  Cartesian 
type  machines.  In  contrast,  the  application  of  a  hexapod  coordinate  measuring  machine 
with  parallel  geometry  is  a  recent  development.  The  calibration  techniques  currently  in 
use  for  finding  the  parametric  errors  for  parallel  structural  machines  require  highly 
accurate  external  measurement  devices  and  the  process  is  very  time  consuming. 
Moreover,  the  calibration  technique  is  different  from  machine  to  machine.  Even  if  the 
kinematic  parameters  have  been  obtained,  these  values  will  keep  changing  when  the 
environmental  temperature  varies.  For  instance,  a  one  meter  long  steel  bar  will  have  a 
length  change  of  13  to  21  urn  (520  to  840  uin)  when  it  is  subjected  to  1  °C  (1.8  °F ) 
temperature  change.  In  the  traditional  Cartesian  machines  one  has  to  perform  a  correction 
procedure  to  compensate  for  the  thermally  induced  errors.  The  development  of  a  fast 
simple,  and  reliable  calibration  method  for  the  HCMM  would  allow  for  the  online  self- 
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calibration  of  kinematic  parameters  and  absolute  measurement  lengths  to  compensate  for 
basic  structural  and  thermally  induced  variations  of  the  machine  elements. 

The  research  presented  in  this  paper  is  dedicated  to  the  development  and  analysis 
of  a  self-calibration  method  for  the  kinematic  parameters  of  the  HCMM.  In  addition,  a 
basic  understanding  of  the  thermally  induced  error  on  the  machine  is  within  the  scope  of 
this  research.  Thus,  the  required  frequency  of  re-calibration  when  the  HCMM  is 
subjected  to  a  thermally  variant  environment  can  be  determined. 


CHAPTER  2 
LITERATURE  REVIEW 

In  general  machine  geometry  error  calibration  includes  four  main  steps:  (1) 
kinematic  modeling  (2)  measurement  (3)  error  identification  (4)  error  compensation. 
Kinematic  modeling  establishes  a  mathematical  model  which  relates  the  machine 
kinematic  parameters  to  the  resulting  pose  of  the  machine.  Measurement  collects 
physical  data  which  contains  information  relating  the  input  of  the  model  to  the  output.  If 
the  geometric  errors  are  directly  measured,  the  kinematic  parameters  of  the  machine  are 
identified  individually.  If  only  the  spatial  coordinates  of  the  end-effector  or  the  probe  tip 
are  measured,  the  geometric  errors  and  kinematic  parameters  have  to  be  identified  by  an 
identification  procedure.  When  the  errors  of  kinematic  parameters  are  identified,  use  of 
this  information  to  correct  the  machine  performance  is  called  error  compensation. 

We  separate  machines  into  two  categories:  serial  kinematic  machines  and  parallel 
kinematic  machines.  The  serial  mechanisms  have  a  open  loop  structure  such  as  machine 
tools,  Cartesian  type  coordinate  measuring  machines,  or  serial  robots.  The  parallel 
mechanisms  have  a  closed  loop  parallel  kinematic  structure  such  as  hexapod  milling 
machines,  Stewart  platforms,  or  the  HCMM.  In  each  category  we  classify  the  calibration 
method  into  two  sub-categories:  "require  external  calibrated  standards"  and  "without 
calibrated  standards"  according  to  the  use  of  the  measurement  device.  The  "without 
calibrated  standards"  calibration  also  is  referred  to  as  the  self-calibration  method.  For 
example,  the  laser  interferometer  measurement  device  (an  external,  calibrated  standard) 
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Figure  2-1  The  kinematic  calibration  categorization 

is  used  for  the  machine  tool  calibration.  In  contrast,  the  center  rod  (an  internal  machine 
component,  not  a  calibrated  standard)  is  used  for  the  HCMM  calibration. 

According  to  the  identification  method,  the  "require  calibrated  standards"  group  in 
the  open  loop  machines  is  further  separated  into  two  groups:  "the  direct  measured"  where 
each  of  the  parametric  errors  is  measured  individually,  and  "the  indirect  measured"  where 
the  parametric  errors  are  identified  simultaneously  by  a  mathematical  method.  The 
calibration  categorization  is  shown  in  Figure  2-1.  A  total  of  five  groups  are  categorized 
by  our  classification. 
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2. 1  Open-Loop  Kinematic  Calibration 

Group  1 :  Direct  Measured 

The  first  group  of  machines  in  calibration  are  the  machine  tool  and  Cartesian  type 
CMM.  In  general,  a  kinematic  error  model  will  be  first  created  to  relate  input  command 
and  the  position  of  the  tool  tip  in  a  machine  tool  or  the  position  of  the  probe  in  a  CMM. 
The  function  can  be  written  as: 

SP  =  [T]8e  (2-1) 
8P  represents  the  position  errors  between  the  actual  and  the  desired  position.  [T]  is  a 
transformation  matrix  that  is  formed  by  the  machine  configuration.  5e  is  the  error 
parameters  needed  to  be  identified  such  as  the  displacement,  straightness,  and  rotation 
error  of  a  machine  carriage.  For  instance,  a  3-axis  machine  tool  or  CMM  has  a  total  of  21 
error  parameters  needed  to  be  identified.  Six  parameters  for  each  axis  and  three 
parameters  for  the  squareness  between  axes. 

In  this  group  the  error  parameters  are  identified  by  direct  measurement.  Normally 
the  laser  interferometers  and  straight  edge  are  used  as  the  external  measurement  device. 
The  main  advantage  in  applying  direct  measurement  methods  is  that  well  established 
methods  are  employed.  The  disadvantages  are:  first,  the  need  to  maintain  many  different 
and  expensive  standards  or  instruments;  second,  when  compared  to  other  methods,  the 
close  attention  to  be  paid  to  sign  conventions  of  error  functions  as  well  as  the  need  for  a 
skilled  operator.  Another  disadvantage  is  that  this  method  is  very  time  consuming;  to 
finish  the  complete  measurement  of  error  parameters  on  a  3-axis  machine  might  take  as 
long  as  a  week.  To  overcome  these  disadvantages  a  new  measurement  instrument  was 
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introduced  by  Ziegert  and  Mize  [5]  to  accomplish  the  parametric  error  measurement  of 
machine  tools  in  hours  rather  in  days.  They  used  a  Laser  Ball  Bar  (LBB)  to  measure  the 
coordinates  of  a  number  of  points  lying  along  the  three  sets  of  mutually  perpendicular 
lines  nominally  aligned  with  the  three  machine  axes,  in  the  working  volume  of  a  machine 
tool.  The  coordinates  of  these  points,  as  determined  by  LBB,  are  used  to  determine  the 
parametric  errors  of  the  three  axes  of  the  machine  tool.  Kulkarni  [6]  compared  the 
traditional  measurement  method  and  the  LBB  method  on  a  White  Sunstrand  Series  20 
milling  machine.  The  LBB  measurements  were  found  to  agree  with  the  conventional 
method. 

The  parametric  errors  of  a  machine  will  vary  under  different  thermal  state.  Thus 
the  parametric  errors  have  to  be  measured  in  different  thermal  states.  Because  of  the 
difficulty  in  linking  the  thermal  state  of  the  machine  to  an  error  function,  only  a  few 
machine  tools  have  taken  full  advantage  of  the  software  compensation  technique. 
Donmenz  [7]  implemented  a  software  compensation  system  on  a  two-axis  turning  center 
to  modify  the  commands  generated  by  the  machine  controller.  In  his  work  the  thermal 
effect  was  also  considered  as  an  error  source. 

In  contrast,  the  software  compensation  technique  has  been  widely  applied  on  the 
CMM.  Since  the  CMM  typically  functions  in  a  tightly  controlled  thermal  environment,  it 
has  less  thermal  variance  than  the  machine  tool.  Therefore,  the  error  function  is  easier  to 
establish. 

Group  2:  Simultaneously  Identified  Using  External  Calibrated  Standards 

For  machines  in  the  second  group  the  error  parameters  are  identified  by  using  a 
best  fit  technique.  For  instance,  the  calibration  of  a  serial  robot  is  in  this  group.  A  typical 
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robot  model  locates  the  link  coordinate  origins  and  defines  transformations  between 
neighboring  links  using  four-parameter  Denavit-Hartenberg  (DH)  representation.  In 
some  cases  the  DH  representation  is  not  sufficient  to  model  the  kinematic  geometry;  for 
instance,  the  DH  representation  is  not  designed  to  represent  a  small  misalignment  in  two 
consecutive  parallel  axes.  A  different  modeling  approach  is  proposed  in  work  by 
Mooring  [8].  This  model  is  based  on  a  set  of  screw  theory  equations  which  is  sometimes 
referred  to  as  Rodrigues'  equations. 

Craig  [9]  demonstrated  the  use  of  the  DH  technique  for  modeling  a  serial  link 
manipulator.  The  procedure  consists  of  assigning  a  coordinate  system  to  each  link  at  the 
joint  axis  and  then  expressing  the  relationship  between  consecutive  coordinate  systems 
with  homogeneous  transformation  matrices.  All  the  individual  link  transformation 
matrices  may  then  be  multiplied  together  to  produce  one  transformation  relating  the 
coordinate  system  in  the  end  effector  to  the  base  coordinate  frame.  This  procedure  is 
called  forward  kinematic  modeling. 

In  contrast,  given  a  desired  end-effector  position  and  orientation  to  find  the 
associated  joint  angle  or  displacement  is  called  inverse  kinematic  modeling.  To  drive  a 
robot  to  a  desired  position,  the  inverse  kinematic  model  has  to  be  used  to  find  the 
corresponding  joint  movement.  If  there  is  any  error  in  the  assumed  kinematic  parameters, 
the  calculated  joint  angles  or  displacements  will  not  drive  the  end-effector  to  the  desired 
position.  Therefore,  an  external  spatial  coordinate  measurement  device  is  usually  used 
during  a  robot  kinematic  calibration  to  measure  the  actual  end-effector  position.  This 
device  may  be  a  precise  CMM  or  a  laser  tracking  system.  A  detailed  description  of  these 
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measurement  devices  is  provided  by  Mooring  et  al.  [10].  When  the  actual  end-effector 
position  is  measured  by  external  devices,  one  can  form  a  residual  equation  as  follows: 
ej=Pactj-f(<lj,k)     7  =  1,2,...,  m  (2-2) 

where  Pact,  j  is  the  measured  end-effector  position  of  the  robot  at  pose  j.     is  the  joint 
readings  at  pose  j.  k  is  the  set  of  kinematic  parameters  needed  to  be  identified. 
Therefore,  a  least-squares  method  can  be  applied  to  estimate  the  value  of  k  that 
minimizes  the  function 

m 

£  =  (2-3) 
;=> 

The  selection  of  the  measurement  configuration  during  robot  calibration,  then, 
plays  an  important  role  in  determining  the  accuracy  and  speed  of  convergence  of  the 
least-squares  identification  algorithms.  The  kinematic  parameter  errors  are  not  equally 
observable.  The  visibility  of  each  unknown  parameter  varies  from  one  robot  to  another. 
There  may  even  be  configurations  in  which  some  of  the  kinematic  parameters  are  not 
observable  at  all.  When  the  unmodeled  error  sources  and  measurement  error  exist,  it 
becomes  impossible  to  obtain  the  exact  value  of  k.  In  minimizing  the  effect  of  the 
inevitable  noise  on  the  results  of  the  estimation,  the  selection  of  robot  measurement 
configurations  may  become  important.  In  other  words,  we  have  to  search  for  the  optimal 
robot  measurement  configurations  which  are  least  sensitive  to  the  measurement  noise  and 
will  obtain  the  best  estimation  results.  Driels  and  Pathre  [11]  investigated  the  factors  that 
will  affect  the  calibration  results.  The  factors  include  number  of  measurements,  the  range 
of  motion  of  the  joints  during  observations,  the  encoder  resolution  and  uncertainty,  and 
the  selection  of  configurations.  In  their  calibration  experiment  simulations  the  condition 
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number  is  used  to  be  an  indicator  of  the  observability  of  the  parameters  to  be  identified. 
Zhuang,  Wang  and  Roth  [12]  adopted  a  simulated  annealing  (SA)  approach  to  obtain 
optimal  or  near  optimal  measurement  configurations  for  robot  calibration.  Zhuang,  Wu 
and  Huang  [13]  developed  another  selection  algorithm  called  "genetic  algorithm"  for 
optimal  measurement  configurations.  In  both  of  the  two  algorithms  the  objective  is  to 
select  the  optimal  configurations  that  minimize  the  condition  number.  Bonn  and  Menq 
[14]  developed  another  index  called  the  observability  index  for  optimal  pose  selection. 
Based  on  this  index,  the  optimal  measurement  configurations  for  kinematic  model 
calibration  are  determined  for  PUMA  type  serial  robots.  Nahvi  and  Hollerbach  [15,  16] 
proposed  two  new  indices:  minimum  singular  value  and  noise  amplification  index  that  are 
used  for  the  optimal  configuration  selections.  They  claim  that  the  noise  amplification 
index  is  more  sensitive  to  calibration  error  than  the  published  observability  indices. 

The  last  step  on  the  robot  calibration  is  to  correct  the  parameter  errors  in  the 
kinematic  model.  Because  the  accuracy  requirement  of  robots  is  much  less  than  the 
machine  tool  or  CMM  and  the  thermal  effect  is  usually  ignored,  robot  error  compensation 
is  much  easier  to  accomplish  than  machine  tool  or  CMM  calibration. 
Group  3:  Self-Calibration 

Because  the  kinematic  calibrations  of  groups  1  and  2  require  an  accurate  external 
measurement  device,  an  experienced  operator,  and  a  large  measurement  time  (between  30 
and  60  hours),  more  and  more  researchers  are  seeking  new  methods  of  calibration  that 
are  less  costly  and  less  time  consuming.  The  self-calibration  technique,  then,  is  starting 
to  be  used  in  some  machines. 
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The  idea  of  the  self-calibration  in  CMM  is  using  a  stable  artifact  such  as  a  ball 
plate  or  ball  array.  During  the  calibration  the  CMM  measure  the  3-D  position  of  balls 
fixed  on  the  plate  for  different  locations  of  the  ball  plate  within  the  measuring  volume  of 
the  CMM.  These  measurements  of  the  ball  plate  in  different  work  space  locations  are 
compared  with  an  off-line  computer  and  kinematic  parameters  are  identified  by  using  the 
least-square  algorithm  that  minimizes  the  measurement  differences  observed  in  various 
locations  throughout  the  work  space. 

The  requirement  for  the  artifact  is  that  it  has  to  be  very  stable.  However,  the 
dimensions  of  the  artifact  need  not  to  be  known  accurately.  The  advantage  is  that  the 
self-calibration  can  easily  be  carried  out  by  non-specialist  CMM  users  with  a  minimum 
setup  time.  Meanwhile,  the  machine  can  easily  be  re-calibrated  after  a  period  of 
operation  or  maintenance  service  and  the  machine  error  in  the  whole  working  volume  can 
be  determined. 

Kruth,  Vanherck  and  Jonge  [17]  performed  the  self-calibration  technique  by  using 
a  ball  plate  on  a  Mitutoyo  FN905  CMM.  The  results  of  the  evaluation  show  that  the 
systematic  geometrical  errors  are  greatly  reduced.  Zhang  and  Zang  [18]  used  a  1-D  ball 
array  to  identify  the  21  parameters  of  a  3-D  CMM.  Belforte  et  al.  [19]  used  a  ball  frame 
to  identify  and  correct  the  parametric  errors  on  nine  CMM  machines.  They  proved  that 
the  self-calibration  technique  can  improve  the  machine  performance  by  a  factor  ranging 
between  three  and  ten.  The  total  hours  required  for  these  calibration  procedures  range 
from  10  to  30  hours  which  is  much  less  than  the  direct  measurement  method.  Pahk  and 
Burdekin  [20]  presented  a  new  parameter  identification  scheme  that  uses  the  stylus  of  the 
CMM  to  measure  a  horizontal  plane  and  obtain  a  locus  of  the  stylus.  They  then  identified 
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the  error  parameters  by  minimizing  the  data  deviation  from  a  mean  plane.  A 
comprehensive  comparison  between  the  direct  method  and  self-calibration  technique  for 
machine  tools  and  CMMs  is  addressed  by  Sartori  and  Zhang  [21]. 

In  contrast,  the  self-calibration  procedure  in  robots  is  much  simpler  than  the 
CMM.  The  idea  of  the  open-loop  robot  self-calibration  is  to  form  a  closed  kinematic 
chain  by  interacting  with  the  environment.  In  other  words,  no  precise  artifact  is  needed  in 
the  robot  self-calibration.  Bennett  and  Hollerbach  [22]  investigate  the  self-calibration  of 
a  6-DOF  manipulator  grasping  a  door  with  a  hinge  joint.  They  show  that  without  any 
external  measurement  device  by  only  using  joint  angle  reading  and  self-motions, 
consistency  conditions  can  be  utilized  to  identify  the  kinematic  parameters.  Khoshzaban 
and  Sassani  [23]  calibrated  a  teleoperated  excavator  with  unsensed  joints  by  adding  an 
additional  linkage  (called  a  calibrator  by  the  authors)  with  some  sensed  joints  to  form  a 
closed  loop.  Zhuang  [24]  shows  that  the  closed  loop  did  not  need  to  involve  physical 
linkages,  but  could  be  formed  by  optical  paths  as  virtual  limbs  to  a  stereo  camera  system. 
He  demonstrated  that  an  uncalibrated  stereo  camera  system  could  be  simultaneously 
calibrated  with  an  uncalibrated  manipulator. 

2.2  Closed-Loop  (Parallel  Mechanism)  Kinematic  Calibration 

The  majority  of  today's  industrial  machines  are  serial  chains  (or  open-loop) 
mechanism.  Cartesian  type  machine  tools  have  been  developed  and  used  for  nearly  200 
years.  The  CMM  and  industrial  robot  calibration  techniques  also  have  many  decades  of 
practical  experience  behind  them.  The  techniques  for  manufacturing  and  calibrating 
these  machines  have  become  very  well  understood  and  mature.  However,  maintaining 
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these  cantilevered  machines  to  be  very  stiff  and  accurate  is  very  difficult  and  expensive. 
In  addition,  most  the  machines  only  provide  single  mode  functionality.  To  make  a  part 
that  requires  milling,  cutting,  welding  and  other  processes  will  require  many  different 
types  of  machines.  In  contrast,  the  parallel  (or  closed  loop)  structures  offer  the  benefits 
of  multiple  functionality  in  a  single  machine,  increased  stiffness  of  the  mechanism,  low 
weight,  and  low  effective  inertia  of  moving  components.  A  performance  comparison 
between  serial  and  parallel  link  manipulators  has  been  addressed  by  El-Khasawneh  and 
Ferreira  [25].  These  advantages  of  the  parallel  type  machines  offer  the  potential  for  better 
performance  than  the  conventional  machines.  As  a  consequence,  parallel  link 
manipulator  concepts  are  increasingly  being  used  in  the  development  of  a  new  generation 
of  advanced  machine  tool. 

A  typical  parallel  type  machine  is  the  hexapod  machine.  The  earliest  known 
hexapod  machine  was  the  tire  tester  designed  by  Gough  in  1949.  However,  this  structure 
is  widely  regarded  as  the  Stewart  Platform  since  Stewart  [26]  popularized  it  through  its 
use  as  a  flight  simulator.  Many  unsuccessful  attempts  were  made  between  1970  and  1990 
to  make  practical  machines.  The  cost  of  computing  strut  lengths  has  until  recently  been 
too  high  for  most  applications.  Now  that  the  cost  of  computing  has  fallen  dramatically, 
many  companies  are  offering  hexapod  based  machines,  such  as  Ingersoll,  Gidding  & 
Lewis,  Hexel,  Geodetics,  ITIA  and  Polytec  PI. 

Group  4:  Simultaneously  Identified  Using  External  Calibrated  Standards 

Since  the  special  demand  for  using  Stewart  platform  as  a  machining  center  was 
only  recently  addressed,  many  of  its  characteristics  are  still  developing.  Although  the 
Stewart  platform  has  many  advantages,  the  disadvantages  include  a  complex  work  space, 
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a  small  orientation  range  and  low  structural  damping.  Because  of  the  complex  structure, 
it  is  necessary  to  analyze  the  contribution  of  each  individual  error  source  to  the  overall 
machine  accuracy.  Before  doing  this,  one  has  to  model  the  forward  and  inverse 
kinematics.  Given  a  set  of  leg  lengths  to  determine  the  platform  position  and  orientation 
is  the  forward  kinematics  problem.  This  problem  could  lead  to  multiple  solutions. 
Merlet  [27]  developed  an  algorithm  to  solve  the  Stewart  platform  forward  kinematic 
problem,  in  which  he  found  that  the  upper  bound  on  number  of  solutions  is  1 320. 
Raghavan  [28]  found  that  a  general  Stewart  Platform  has  40  real  solutions.  However, 
Wang  and  Masory  [29]  presented  a  numerical  method  for  solving  the  forward  and  inverse 
kinematics  and  analyzed  the  effect  of  the  manufacturing  tolerance,  installation  errors  and 
link  offsets  on  the  platform  accuracy.  Wang,  Masory  and  Zhuang  [30]  also  established  a 
calibration  algorithm  based  on  the  error  model  derived  by  forward  kinematics. 

Unlike  the  forward  kinematics  problem,  the  inverse  kinematics  problem  is  easy  to 
solve,  and  is  defined  as,  "Given  the  platform  position  and  orientation,  compute  the 
manipulator  leg  lengths."  Because  of  the  difficulty  in  the  forward  kinematic  modeling, 
many  researchers  used  the  inverse  kinematics  to  find  the  error  sources  and  how  the  error 
sources  affect  on  the  accuracy  of  the  parallel  machines.  According  to  the  analysis  of 
Ropponen  and  Arai  [31]  and  Patel  and  Ehmann  [32],  the  error  model  of  a  Stewart 
platform  based  machine  may  be  expressed  as 

8A  =  jm  +  NSA  (2-4) 
where  5A  is  the  leg  length  errors.  811  is  the  platform  position  and  orientation 
errors.  6A  is  the  joint  position  errors.  J  is  the  Jacobian  matrix  of  the  Stewart  platform.  N 
is  a  transformation  matrix  mapping  the  joint  error  to  leg  length  error.  Patel  and  Ehmann 
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completed  an  automated  error  analysis  software  called  SPAEA  (Stewart  Platform 
Automated  Error  Analysis)  for  platform  accuracy  analysis.  Assuming  the  error 
components  are  given,  a  platform  position  and  orientation  accuracy  analysis  can  be 
performed  using  the  nominal  parameters  in  the  SPAEA.  The  idea  is  to  invert  Equation  2- 
4  and  obtain  the  platform  error  function 

m  =  J-l(SA-N8A)  (2-5) 

The  same  error  model  as  shown  in  Equation  2-4  also  has  been  developed  by 
Soons  [33]  about  the  same  time  as  Patel  and  Ehmann.  He  also  used  the  inverse 
kinematics  idea  but  he  derived  the  error  functions  more  concisely  than  Patel  and  Ehmann. 

Another  calibration  technique  was  suggested  by  Zhuang  and  Roth  [34]  who 
proposed  an  idea:  by  fixing  one  leg  at  a  time  and  moving  the  other  five  legs  during 
measurement  process,  one  is  able  to  compute  kinematic  parameters  one  leg  at  a  time. 
This  approach  requires  relatively  little  computation.  However,  the  identified  parameters 
may  not  be  optimal  because:  (1)  The  measurable  workspace  is  small  because  of  the 
restricted  motion  pattern  of  the  manipulator;  (2)  Kinematic  parameters  of  individual  links 
are  computed  separately,  therefore,  the  coupling  effects  among  the  legs  are  not  fully 
explored;  and  (3)  The  model  used  did  not  include  all  the  kinematic  parameters  of  the 
platform  because  they  assumed  the  joints  were  perfect. 

Once  the  error  model  has  been  completed,  the  next  step  for  kinematic  calibration 
is  the  measurement.  An  external  measurement  device  has  to  be  used  for  performing  this 
task.  This  device  measures  the  position  and  orientation  errors  of  the  platform  shown  in 
Equation  2-5.  Soons  [35]  developed  a  computer-aided  experimental  design  technique  for 
the  configuration  selection  that  improves  the  condition  of  error  parameter  identification. 
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Since  the  observed  errors  are  the  linear  combination  of  the  error  components  that  are 
shown  in  Equation  2-5,  the  least-squares  method  is  used  to  identify  the  error  parameters 
simultaneously. 

In  addition  to  the  hexapod  machine,  Maurine  and  Dombre  [36]  presented  a 
calibration  procedure  for  a  4-DOF  parallel  robot.  They  used  an  inexpensive  laser 
displacement  sensor  installed  on  the  robot  platform  to  measure  an  external  grid  pattern 
and  several  cylinders  placed  in  a  circle  on  a  work  table.  They  compared  the  errors 
between  the  values  calculated  from  kinematic  model  and  the  values  from  measurement. 
The  least-squares  method  then  was  used  for  identifying  the  error  parameters.  Wallack  and 
Mazon  [37]  used  a  calibration  peg  and  multiple  crossbeam  sensors  on  a  4-DOF  robot  to 
identify  the  rotation  axis  and  translation  axis. 
Group  5:  Self-Calibration 

As  we  mentioned  earlier,  when  the  kinematic  calibration  needs  an  external 
measurement  device,  this  process  will  require  experienced  operators  and  be  very  time 
consuming.  Thus,  the  self-calibration  techniques  are  an  attractive  alternative.  This 
technique  just  needs  some  redundant  sensors  or  kinematic  constraints  and  moving  the 
machine  to  different  configurations,  then  the  error  parameters  can  be  identified.  Though 
the  self-calibration  can  be  simply  performed,  only  a  few  papers  have  been  published  at 
this  time. 

Zhuang  and  Liu  [38]  installed  several  redundant  sensors  on  the  U-joints  of  a 
Stewart  platform  based  machine,  and  identified  the  30  error  parameters.  The  idea  is  to 
adjust  the  nominal  values  to  minimize  the  difference  between  the  measured  values  from 
redundant  sensors  and  the  calculated  values  from  the  forward  and  inverse  kinematic 
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models.  Nahvi,  Hollerbach  and  Hay  ward  [15]  presented  an  autonomous  kinematic 
calibration  method  for  a  3-DOF  shoulder  joint  parallel  robot.  They  used  the  features  of 
the  redundant  kinematics  that  form  multiple  closed  loops,  thus,  adding  up  all  vector 
components  in  each  closed  loop  should  yield  a  zero  vector.  If  the  nominal  parameters  are 
incorrect,  however,  the  closed  loop  vector  summation  will  yield  a  residual  error  vector. 
The  kinematic  parameters,  then,  are  identified  by  minimizing  the  residual  vectors.  They 
applied  the  smallest  singular  value  along  with  the  condition  number  of  the  identification 
Jacobian  as  the  pose  selection  strategy  that  yields  a  maximum  observability  of  the 
parameters.  Hollerbach  and  Lokhorst  [39]  used  the  same  idea  in  [15]  and  presented 
experimental  results  on  the  parameter  identification  of  a  RSI  6-DOF  hand  controller,  a 
two-loop  parallel  mechanism  comprised  of  three  6-DOF  arms  with  potentiometers  on  the 
first  joint  of  each  arm. 

2.3  Conclusion 

Based  on  the  above  review,  if  the  self-calibration  technique  is  not  adopted,  then 
expensive  special  instrumentation  is  required  to  perform  the  parameter  identification  for 
the  HCMM.  Even  if  the  instrumentation  is  available,  specially  trained  technicians  are 
needed  to  make  the  measurements.  Large  data  collection  times  are  required  to  measure 
the  parametric  errors  of  machines.  The  cost  and  difficulty  of  the  external-artifact  based 
method  prevents  them  from  being  used  to  re-calibrate  frequently.  In  contrast,  the  self- 
calibration  technique  is  less  expensive,  faster,  and  easier  to  operate.  It  allows  the 
machine  to  be  re-calibrated  more  frequently.  Hence,  it  would  be  a  better  alternative  to 
identify  the  parametric  errors  of  the  HCMM  by  a  self-calibration  technique. 


CHAPTER  3 

MODELING,  MEASUREMENT  AND  IDENTIFICATION 

As  mentioned  earlier,  the  laser  interferometric  displacement  measurement  system 
is  installed  in  each  of  the  six  struts  and  measures  the  length  changes  of  each  strut. 
However,  in  order  to  calculate  the  coordinates  of  the  probe  of  the  HCMM,  it  is  necessary 
to  know  the  absolute  length  of  each  strut  and  the  distance  between  the  three  base  spheres. 
The  nominal  lengths  could  be  obtained  by  using  a  conventional  linear  measurement 
device  such  as  a  scale.  However,  to  find  the  kinematic  parameter  values  with  sufficient 
accuracy,  a  calibration  procedure  is  necessary. 

The  process  of  kinematic  calibration  consists  of  the  modeling,  measurement,  and 
identification  steps.  Modeling  refers  to  the  choice  of  a  functional  relationship  between 
the  kinematic  parameters  and  the  resulting  poses  of  the  center  rod  of  the  HCMM.  The 
model  selected  should  account  for  all  the  factors  considered  to  be  significant  in 
contributing  to  HCMM  accuracy. 

Physical  data  are  then  collected  from  moving  the  center  rod  into  different  poses. 
This  data  set  contains  information  relating  the  input  of  the  model  (the  reading  of  strut 
displacements)  to  the  output  of  the  model  (the  center  rod  length). 

The  mathematical  process  of  using  the  data  collected  to  identify  the  parameters  of 
the  model  is  the  third  step  in  calibration.  It  is  important  to  make  the  identification 
process  fast  and  accurate.  In  this  chapter  we  will  describe  these  three  calibration 
processes. 
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3.1  Kinematic  Modeling 

The  first  step  in  any  calibration  procedure  is  to  obtain  a  valid  kinematic  model. 
The  purpose  is  to  relate  the  outputs  of  joint  displacement  to  the  manipulator  pose.  There 
are  two  basic  forms  that  any  kinematic  model  may  take.  The  forward  model  computes 
the  manipulator  pose  from  the  given  joint  displacement  readings.  The  inverse  model,  on 
the  other  hand,  is  used  to  find  the  set  of  joint  displacements  when  a  manipulator  pose  is 
given.  When  we  formulate  the  kinematic  model  of  the  HCMM,  the  following 
assumptions  are  made: 

1 .  The  spherical  joints  produce  perfect  spherical  motion  and  the  centers  of  the  reference 
spheres  are  their  exact  rotation  centers. 

2.  No  noticeable  thermal  effect  occurs. 
The  World  Coordinate  System  Assignment 

We  define  a  fixed  world  coordinate  system,  to  which  all  coordinate-dependent 
kinematic  parameters  are  referred.  In  general  it  can  be  placed  anywhere  in  space.  A 
principle  for  placing  this  coordinate  system  is  not  to  introduce  any  redundant  kinematic 
parameters.  In  this  case  the  origin  of  the  coordinate  system  is  placed  on  one  of  the  base 
sphere  centers,  BS1,  as  shown  in  Figure  3-1. 

The  x  axis  lies  along  a  line  connecting  two  base  sphere  centers.  The  plane  passing 
through  the  centers  of  three  base  spheres  is  chosen  as  x-y  plane.  The  z  axis  is  along  the 
normal  direction  of  the  x-y  plane.  The  BS2  has  coordinate  of  (r,0,0)  and  BS3  is  located  at 
(b,h,0).  Therefore,  in  this  coordinate  system  setup,  a  total  of  nine  parameters  are  required 
to  determine  the  coordinates  of  the  upper  and  lower  sphere  centers,  Pup  and  Piow. 
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P_up 


P_low 


Figure  3-1  The  HCMM  coordinate  system 

Forward  Kinematic  Model 

The  constant  center  rod  length  provides  the  kinematic  chain  which  constrains  the 
motion  of  the  apexes  of  the  two  tetrahedrons.  Its  length  is  defined  as  the  distance 
between  the  two  sphere  centers.  If  the  locations  of  the  three  base  spheres  and  the  absolute 
length  of  each  strut  are  given,  the  unique  position  and  orientation  of  the  center  rod  can  be 
determined. 

Because  of  the  telescoping  strut  design,  each  of  the  six  struts  may  change  its 
length  to  satisfy  the  kinematic  constraint  when  the  center  rod  is  moved  from  one  position 
to  the  other.  We  define  the  "length  change"  as  the  length  difference  between  the  initial 
and  the  current  configuration.  For  the  ith  configuration,  the  strut  length  is  the  sum  of  the 
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initial  (or  starting)  length  and  its  length  change  at  the  current  configuration.  Let's  say  the 

kth  strut,  with  its  initial  length  of  a(k),  has  a  length  change  of  S(k,  i)  at  the  ith 

configuration.  At  this  instance  its  absolute  length  is  expressed  as: 

l(k,i)  =  a(k)  +  S(k,i) 
where 

t-U  .6  (3"1) 

i  —  l,..,m 

From  the  above  expression  it  is  noted  that  the  initial  length  of  a(k)  is  a  constant 
value  and  the  length  change  of  S(k,  i)  will  vary  from  configuration  to  configuration.  The 
centers  of  the  upper  and  lower  spheres,  Pup  and  Piow,  can  be  obtained  from  Equations  3-2 
to  3-7. 


X.(i)  = 


/2(l,0-/2(3,Q  +  r2 
2r 


(3-2) 


YUH)  = 


1 2 (3,  Q  -  I2 (5,  Q  +  2XU (i)[r  -b]-[r2-b2  -h2] 

2h 


(3-3) 


(3-4) 


X,(i)  = 


l2(2,i)-l2(4J)  +  r2 
2r 


(3-5) 


Y,(i)  = 


1 2  (4,  Q  -  I2  (6,  Q  +  2X,  (i)[r  -b]-[r2  -b2  -h2] 

2h 


(3-6) 


(3-7) 


Therefore,  the  center  rod  length  can  be  calculated  as: 


Lc=V[*,(0-X,(i)]2  +[Yu(i)-Yl(i)f  +[Zu(i)-Zl(i)f 


(3-8) 
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Inverse  Kinematic  Model 

In  the  forward  kinematic  model  the  strut  displacements  are  treated  as  inputs.  The 
HCMM  controller,  on  the  other  hand,  may  receive  a  desired  center  rod  pose  as  the  input 
and  must  compute  the  joint  displacements  necessary  to  achieve  this  pose.  In  other  words, 
the  inverse  kinematic  model  is  required  for  use  in  the  controller. 

Because  of  the  simple  geometry  of  our  machine,  the  inverse  model  can  be  easily 
obtained  as  the  following, 


lxfi+Y;+Z;-a(\) 

(3-9) 

S(2)  =  . 

Jxf  +  Y,2+Zl  -a(2) 

(3-10) 

S(3)  =  , 

l{Xu-rf  +  Y*+Zl-aO) 

(3-11) 

S(4)  =  . 

J(X,-r)2+Y?+Zf  -a(4) 

(3-12) 

S(5)  =  , 

](Xu-b)2+(Ytt-h)2+Z2u  -a(5) 

(3-13) 

5(6)  =  1 

/(X, -6)2+(y;-/i)2+Z,2-a(6) 

(3-14) 

3.2  Measurement 


The  fixed  length  center  rod  and  the  closed  loop  kinematics  of  the  HCMM  provide 
a  very  useful  feature  in  the  calibration  procedure.  Unlike  the  conventional  calibration 
technique,  one  can  use  the  kinematic  constraints  and  the  displacement  reading  of  each  of 
the  struts  to  perform  the  kinematic  calibration.  Therefore,  this  calibration  procedure  is 
called  "self-calibration"  since  no  external  measurement  device  is  required. 
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For  the  calibration  of  a  conventional  CMM,  the  measurement  process  is  time 
consuming,  laborious,  and  prone  to  human  error.  In  contrast,  the  measurement  process 
for  HCMM  is  simply  moving  the  center  rod  to  a  number  of  different  positions  within  its 
work  volume,  and  recording  the  displacement  reading  in  each  strut.  Thus  one  can  collect 
data  very  easily.  However,  it  is  necessary  to  minimize  the  number  of  measurement, 
because  we  assume  that  there  is  no  thermal  effect  during  the  calibration  process.  Taking 
too  many  measurements  could  create  more  heat  released  by  the  actuator  and  cause 
significant  thermal  deformation  in  the  structure. 

We  will  show  that  the  self-calibration  process  is  very  sensitive  to  the  level  of 
measurement  noise.  The  measurement  noise  may  include  the  modeling  imperfection  or 
measurement  error  from  the  laser  interferometer.  To  eliminate  or  minimize  the  effect 
owing  to  measurement  noise,  then,  becomes  one  of  the  most  difficult  problems  in  the 
parameter  identification.  In  fact,  it  is  impossible  to  include  every  possible  error  source  in 
a  kinematic  model  or  to  completely  eliminate  the  measurement  error.  When  the  un- 
modeled  error  sources  and  measurement  error  exist,  it  becomes  impossible  to  obtain  the 
exact  values  of  the  parameters  used  in  the  model.  In  general,  the  magnitude  of  the  noise 
itself  is  quite  small,  but  its  effect  on  the  estimation  accuracy  may  be  rather  significant.  To 
minimize  the  effects  of  the  inevitable  noise  on  the  results  of  estimation,  the  selection  of 
measurement  configurations  may  become  important.  The  selection  of  the  optimal 
measurement  configuration  for  calibration  will  be  discussed  in  a  later  chapter. 
Error  Sources 

As  we  have  mentioned  the  self-calibration  method  is  very  sensitive  to  the 
measurement  noise  (or  error),  it  is  necessary  to  realize  and  analyze  the  error  sources  and 
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to  reduce  their  magnitude  during  the  design  and  manufacturing  stages.  The  error  sources 
could  be  (1)  the  displacement  measurement  error  of  the  laser  interferometer  (2)  the 
imperfections  of  the  spherical  joint  and  the  spheres  (3)  the  thermal  deformation  of  the 
structure.  A  detailed  description  on  the  error  sources  will  be  discussed  in  Chapter  7. 


The  mathematical  process  of  using  the  data  collected  to  identify  the  parameters  of 
the  model  is  the  third  step  in  the  calibration.  Because  of  the  incorrect  nominal  values, 
measurement  noise  and  modeling  error,  the  calculated  center  rod  length  from  the  model 
would  not  have  the  same  value  at  each  pose.  If  the  true  center  rod  length,  Lctrue ,  could  be 
obtained,  the  residual  equation  of  the  center  rod  could  be  presented  by  Equation  3-15. 


where  Lc,  is  the  predicted  center  rod  length  calculated  by  the  model  at  the  ith 
configuration.  Thus  the  kinematic  parameters  can  be  obtained  by  adjusting  their  values  to 
minimize  the  residual  of  the  center  rod  length.  However,  in  order  to  implement  this 
scheme  it  is  necessary  to  know  the  true  center  rod  length,  or  at  least  the  length  that  is 
precise  enough  for  the  calibration  purposes.  It  will  be  a  difficult  problem  to  find  a 
measurement  device  and  technique  that  can  be  performed  quickly  and  easily  to  give  us  a 
sufficiently  accurate  value  for  the  center  rod  length. 

Another  identification  algorithm  to  substitute  Equation  3-15  can  be  written  as 
follows: 


3.3  Identification 


SLCi  =  Lc,rue  ~ 


Lc 


(3-15) 


5Lcu  =  Lc,.  -  Lc  j 


(3-16) 
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A  A 

The  Lc,  and  Lc ,  are  the  calculated  center  rod  lengths  in  ith  and  jth  configurations 

respectively.  We  call  the  dLc^  the  relative  residual  of  the  center  rod  length.  Using 
Equations  3-1  to  3-8,  the  relative  residual  may  be  rewritten  as 

SLctj=  f(a,r,b,h,St)- f(a,r,b,h,Sj)  (3-17) 

where  Si  is  a  set  of  strut  displacements  associated  with  configuration  i.  Equation  3-17  is 
a  feasible  algorithm,  but  to  perform  the  mathematical  identification  would  be  very 
computationally  intensive.  For  instance  by  taking  m  configurations,  will  generate 

fflX(ffl-l) 

  relative  residual  equations  which  is  1225  for  m  =  50.  To  solve  such  large 

residual  equation  would  require  long  computation  times.  Therefore,  a  modified 
algorithm  called  reference  residual  is  used  in  the  HCMM  calibration  procedure.  It  is 
represented  as  follows: 

A 

ei  =  Lc;  -  Lei2  (3-18) 
where  Lc,  is  an  assigned  value  which  is  an  approximate  guess.  It  is  not  necessary  for  this 
value  to  be  known  exactly  and  it  will  be  treated  as  the  tenth  kinematic  parameter  which 
needs  to  be  identified.  The  objective  is  then  to  find  a  set  of  parameters  including  the  Lct 
to  minimize  the  reference  residual,  ex. 
Least-Squares  Method 

Our  goal  is  to  find  the  optimal  parameters  that  minimize  the  reference  residual. 
There  are  many  optimization  techniques  one  can  use.  One  technique  is  to  find  a  set  of 
parameters  that  minimize  the  sum  of  the  absolute  reference  residual  which  is  expressed  in 
Equation  3-18. 
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Minimize  ^Je,  | 


(3-19) 


i=l 


A  gradient  search  may  be  used  for  finding  the  unknown  parameters.  But  this 
method  is  not  popular  because  it  is  computationally  demanding  and  mathematically 
intractable.  On  the  other  hand,  from  our  literature  review,  the  least-squares  method  was 
seen  as  the  most  popular.  It  has  two  big  advantages.  First,  large  errors  are  heavily 
penalized.  The  other  advantage  is  mathematical  tractability.  The  formula  giving  the 
least-squares  estimates  is  obtained  by  quite  simple  matrix  algebra,  and  the  estimates  are 
computed  as  the  solution  to  a  set  of  linear  equations.  Moreover,  the  properties  of  the 
estimates  are  relatively  easy  to  analyze.  The  least-squares  method  can  be  expressed  as 


Minimize  ^,e,2 


(3-20) 


1=1 


Linearized  Least-Squares  Method  (Gauss-Newton's  Method) 


Since  Equation  3-18  is  highly  non-linear,  one  method  to  solve  it  is  to  linearize  this 
equation.  At  the  ith  configuration  of  HCMM,  by  using  the  Taylor  expansion,  Equation  3- 
21  can  be  written  as 


e;  =e,  .  + 


"4.  da, 


de: 


Aa, +•••+• 

a  a, 

K  6 


Aa6+— J- 
dr 


de: 


db 


Ab 


+  ■ 


de: 


dh 


Ah  + 


de.. 


die. 


(3-21) 


ALcr  +  higher  order  terms 


where  vector  A0  is  the  nominal  value  for  parameters  {a, ,  •  •  • ,  a6 ,  r,  b,  h,  Lcr }  7  . 

Suppose  that  the  center  rod  is  placed  into  m  poses,  this  provides  a  vector  of  e 


which  is  written  as 
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V 

V 

+ 

<9e,     <9e,     <9e,  <9e, 


dh  dLcr 


Aa6 
Ar 
Ah 
Ah 

KALcrJ 


or 


+ higher  order  terms 


E  —  Eo  +  JAA  +  higher  order  terms 


where  J  = 


da, 


Br     dh     dh  dLc 


da, 


dh     dh  dLc, 


AA  = 


(  Aa^ 


Aa6 
Ar 
Ab 
Ah 
KALcrJ 


(3-22) 


(3-23) 


The  values  of  E0  and  J  are  evaluated  at  A0.  The  J  matrix  is  referred  as  the  Jacobian 
matrix.  If  the  nominal  values  of  A0  are  close  enough  to  the  true  values,  the  higher  order 
terms  in  Equation  3-23  then  can  be  ignored.  In  this  case  the  reference  residual  can  be 
written  as 

E  =  E„  +JAA  (3-24) 
Therefore,  the  objective  function  of  Equation  3-20  can  be  rewritten  as 

Minimize  f  ={E0  +  JAA)T(E0  +  JAA)  (3-25) 


35 

This  performance  criterion  is  quadratic  in  the  elements  of  AA.  The  minimum  value  of /is 
obtained  when 

df 

J^  =  °  (3-26) 

together  with  the  condition  that  at  the  value  of  AA  that  satisfies  Equation  3-26,  the 
Hessian  matrix  of /is  positive  semidefinite: 

d2f 

To  differentiate /  the  following  matrix  derivation  formulas  are  used: 

S-(aTHe)  =  He  (3-28) 
a  a 

j^(eTHa)  =  HTe  (3-29) 

—  (aTHa)  =  (H  +  HT)a  (3-30) 
where  H  is  any  square  matrix  and  a  and  e  are  vectors. 

Differentiating  Equation  3-25  and  substituting  into  the  condition  given  in  Equation  3-26 
yields 

df 
dAA 

which  yields 


=  -2JTEo+(JTJ  +  JTJ)AA  =  0  (3-31) 


AA  =  (JTjylJTE0  (3-32) 
During  each  iteration  the  guess  parameter  is  updated  as 

=  Ku  +  AA  (3.33) 
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This  procedure  is  iterated  until  the  largest  value  in  vector  A  is  less  than  the  desired 
tolerance.  This  method  also  is  referred  to  as  the  "Gauss-Newton"  method  or  the  iterative 
linear  least-squares  method.  The  flow  chart  of  linearized  least-squares  method  is  shown 
in  Figure  3-2. 
Gradient  Method 

Without  linearizing  the  non-linear  least-square  equation,  there  are  still  several 
optimization  methods  which  can  be  used.  One  may  try  the  gradient  method.  The 
iterative  technique  is  used  for  finding  the  optimal  parameters.  Let  ay  denote  the  estimate 
of  the  optimal  solution  at  the  beginning  of  y'th  iteration.  The  jth  iteration  consists  of  the 
computation  of  a  search  vector  p_j  from  which  the  new  estimate      is  obtained  according 
to 


If  the  value  of  p_j  equals  to  the  negative  of  the  gradient  of  the  objective  inj'th 
iteration,  this  is  called  the  steepest  descent  direction.  In  general  the  choice  of  vector  p_j  is 
not  necessary  to  along  the  direction  of  steepest  decent.  According  Newton's  method  the 
search  vector  p_s  is  chosen  by 


where  Hs  and  gj  are  the  Hessian  matrix  and  gradient  of  the  objective  function  in  jth 
iteration.  They  can  be  expressed  in  matrix  form  as 


(3-34) 


P.  =  Hrlg 
—j      '  —i 


(3-35) 


8  =~2/V 
—j  ! 


(3-36) 


Hj=2JjTJj+2Qj 


(3-37) 
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Figure  3-2  Flow  chart  for  linearized  least-squares  method 


38 


where  J  is  the  Jacobian  matrix  in  yth  iteration,  gj  is  the  objective  evaluated  at  a  =  a  .j .  Qj 


is  the  value  of 


Qj  =  X£,-(V2[£]) 


(3-38) 


Substituting  Equation  3-36  and  3-37  into  Equation  3-35  we  obtain 


Pj=(JjTJj+QjTiJjTej 


(3-39) 


In  Equation  3-34  if  the  estimate  parameter  is  close  to  the  true  values,  the  Qj  is 
small  enough  to  be  ignored.  Also  assuming  the  Otj  =1,  then  the  identification  algorithm 
reduces  to  the  linearized  least-squares  method  or  the  Gauss-Newton  method. 
Levenberg-Marquardt  (LM)  Method 

The  Gauss-Newton  algorithm  breaks  down  at  the  singularities  of  the  identification 
Jacobian  and  converges  very  slowly  near  such  singularities.  Therefore,  the  LM  method  is 
introduced  as  an  improvement  to  the  Gauss-Newton  method.  The  technique  is  designed 
to  overcome  problems  related  to  singularity  of  the  matrix  fj.  The  idea  is  to  modify 
Equation  3-39  by  adding  a  time  varying  nonnegative  scalar  coefficient,  |ij ,  as  follows: 


where  /  is  an  identity  matrix. 

The  value  jjj  is  initially  set  to  some  positive  value  (say,  fXj  =  0.01).  At  the 
beginning  of  each  iteration,     is  reduced  by  a  factor  of  k  (say,  it  =  10)  in  attempt  to  push 
the  algorithm  closer  to  the  Gauss-Newton  method.  If,  for  the  chosen  value  of  (ij,  no 
reduction  in  the  objective  function  is  achieved,  the  value  of  jij  is  increased  by  the  factor  k. 
For  more  detailed  description  for  selection  of  jLLj  and  k  is  provided  by  Scales  [40]. 


(3-40) 
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3.4  Conclusion 

In  this  chapter  we  developed  a  self-calibration  algorithm  that  is  feasible  for  the 
identification  of  the  HCMM  parameters.  The  idea  for  the  self-calibration  is  to  introduce  a 
kinematic  constraint  so  that  the  HCMM  can  be  calibrated  without  adopting  any  external 
measurement  device.  The  center  rod  length  can  be  treated  as  a  new  parameter  and  its 
value  is  not  necessary  known.  The  optimal  kinematic  parameter  values  can  be  obtained 
by  applying  the  least-squares  method  to  reduce  the  reference  residual  of  the  center  rod. 


CHAPTER  4 
MEASUREMENT  STRATEGY 

The  most  notable  feature  of  a  self-calibration  technique  is  that  no  external 
measurement  device  is  needed.  By  moving  the  center  rod  into  a  number  of  poses  and 
recording  the  displacement  data  of  each  strut,  one  can  determine  the  kinematic  parameters 
of  the  mechanism.  However,  this  self-calibration  technique  is  very  sensitive  to  the  level 
of  measurement  noise.  If  the  measurement  noise  is  too  big,  the  calibration  results  may  be 
very  poor. 

The  measurement  noise  may  include  the  modeling  imperfection  or  measurement 
error  from  the  laser  interferometer.  To  eliminate  or  minimize  the  effect  owing  to 
measurement  noise,  then,  becomes  one  of  the  most  difficult  problems  in  the  parameter 
identification.  In  fact,  it  is  impossible  to  include  every  possible  error  source  in  a 
kinematic  model  or  to  completely  eliminate  the  measurement  error.  When  the  un- 
modeled  error  sources  and  measurement  exist,  it  becomes  impossible  to  obtain  the  exact 
values  of  the  parameters  used  in  the  model.  In  general,  the  magnitude  of  the  noise  itself 
is  quite  small,  but  its  effect  on  the  estimation  accuracy  may  be  rather  significant.  For 
example,  in  some  HCMM  configuration,  center  rod  length  errors  are  so  insensitive  to  the 
error  parameters  that  the  unmodeled  error  sources  become  the  main  causes  of  the  center 
rod  length  error.  When  the  center  rod  length  error  calculated  from  these  configurations, 
the  estimation  result  will  obviously  be  poor.  We  want  the  effect  of  measurement  error 
and  unmodeled  errors  to  be  less  significant  to  the  determination  of  error  parameters, 
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consequently  a  better  estimation  of  parameters  will  be  expected.  Therefore,  to  minimize 
the  effects  of  the  inevitable  noise  on  the  results  of  estimation,  the  selection  of 
measurement  configurations  becomes  very  important. 

Many  researchers  have  proposed  a  variety  of  observability  indices  to  quantify  the 
goodness  of  pose  selection.  These  indices  are  based  on  the  singular  value  decomposition 
(SVD)  of  Jacobian  matrix  of  the  differential  kinematics.  These  indices  are  (1)  condition 
number  (2)  observation  index  (3)  smallest  singular  value  (4)  noise  amplification  index. 
Since  the  condition  number  is  the  most  widely  used  index  for  configuration  selection,  we 
choose  the  condition  number  as  the  index  to  quantify  the  goodness  of  pose  selection.  The 
following  section  will  show  how  the  condition  number  is  derived.  For  more  information 
on  the  other  observability  indices,  one  can  obtain  it  from  Appendix  A. 


Before  starting  any  description  on  the  condition  number,  it  is  better  to  introduce 
the  matrix  norms  and  the  SVD.  An  induced  matrix  norm  for  square  matrix  A  is  defined 
as 


Condition  Number 


AX 


\\X 


for  X  *  0 


(4-1) 


V 


which  satisfies  all  the  following  conditions  that  define  the  norm  in  general 

(1)N|>0,  |H|  =  0  if  only  if  A  =  0 


(4-2) 


(2)M||  =  H||/1 


(4-3) 


(3)||A  +  5||<||A||  +  ||fi|| 


(4-4) 


(4)||Afi||<||A||||fi| 


(4-5) 
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Depending  on  which  vector  norm  is  used  we  will  see  different  norm  values.  For  a 
n-dimension  vectors  with  real  elements  the  most  common  norms  are  1-norm,  2-norm, 


and  infinity-norm  as  follows: 


=Xhl 


(4-6) 


AIL  = 


f  n  \/2 

V  i=i  J 


\\x\\m  -  max  pcf.  I  for  \<i<n 


(4-7) 


(4-8) 


One  can  show  the  1  -norm  produces  the  maximum  column  norm  for  the  induced 
matrix  norm  in  Equation  4-1.  2-norm  gives  the  largest  singular  value  of  A  which  also  is 
the  largest  eigenvalue  of  AT A  .  The  infinity  norm  yields  the  maximum  row  sum.  The  1- 
norm  and  °o-norm  can  be  written  as  follows: 


max  |X  k|    f°r  1  *  J  *  n\ 


(4-9) 


=  max  \  ^ \ajj  for  1  <  /  < , 


(4-10) 


7=1 


Next  the  SVD  will  be  defined.  For  a  mx  n  matrix  A  with  m>  n  (this  assumption 
is  made  for  convenience  only;  all  the  results  will  also  hold  if  m<n),  there  is  a  singular 
value  decomposition  as 

A  =  UTQT  (4-11) 


where 


Z  = 


v0  0, 


(4-12) 
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U  is  a  m  x  m  matrix,  Zi  is  rxr  diagonal  matrix  whose  diagonal  entities  are  non- 
zero singular  values  cr,  >  o~2  >  ■■  >  or  >  0 ,  and  Q  is  nxn  matrix.  The  U  and  Q  are 
orthonormal  matrices.  The  columns  in  U  are  the  left  singular  vectors  of  matrix  A  and  the 
columns  in  Q  are  the  right  singular  vectors  of  A.  If  A  is  nonsingular  or  full  rank  matrix, 
then  r  equals  to  //.  In  this  case  £  can  be  expressed  as 


tr,  0  0 

0  a,  0 

0  o  '•. 

0  0  ••■ 

0  0  0 


0^ 

0 

0 

0 


(4-13) 


v  0     0     0    •••    0  , 

'  mxn 


More  detail  about  how  to  derive  the  SVD  can  be  obtained  form  Klema  and  Laub 

[40]  and  Leon  [41]. 

After  the  matrix  norm  and  SVD  have  been  introduced,  we  start  to  analyze  the 
numerical  sensitivity  of  the  least-squares  identification  problem.  Let  A  be  a  nonsingular 
matrix  of  order  n  and  consider  the  system  Ax  =  b .  The  vector  b  corresponds  to  the 
measured  data  and  as  such  is  subject  to  uncertainty.  The  matrix  A  is  related  to  the 
measurement  configuration  and  the  system  model  and  also  is  subject  to  errors.  Thus: 

x  +  8x  =  (A  +  5AYl(b  +  8b)  (4-14) 
where  the  solution  error  5x  reflects  the  sensitivity  of  the  identification  algorithm  to  errors 
in  the  model  and  the  data  characterized  in  Equation  4-14  as  8A  and  8b. 
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Assume  first  that  8A  =  0  but  8b  *  0.  Then 

IN<IM||<5/>||  (4-15) 

Similar  to  Equation  4-15  one  may  write 

fcNlMHWI  (4-16) 
Combing  Equations  4- 1 5  and  4- 1 6  and  assume  that  b  *  0,  we  find  that 

N    ii  nil    ,11 \M 

^^MIIlA-'l^  (4-17) 

Also  from  Equation  4-14,  assume  that  8A  =  0.  We  have 

8b  =  A8x  (4-18) 
We  may  write  the  inequality  equation  as 

||5fc||<  Milfoil  (4-19) 
From  the  system  equation  Ax  =  b,  we  have 

bMK'liy  (4-20) 
Combining  Equation  4-19  and  4-20,  we  find  that 

i    Ml  INI 


INIIK'll  H  ~  llil 

From  Equation  4- 1 7  and  4-2 1 ,  we  may  write  that 

I      \\5b_l  \\Sx\\ 


<jf<\\A\\\\A-'\\\^  (4-22) 


(4-21) 


IHHK'II  Ik, 

The  number  |U|||a~'  ||  is  called  the  condition  number  of  A  and  will  be  denoted  by 
cond(A).  Thus: 
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1      \\Sb\\  \\8x\ 


(4-23) 


cond(A)  \\b\\  \\x 


The  inequality  Equation  4-23  related  the  size  of  relative  error  t-tt  to  the  relative  residual 


.  If  the  condition  number  is  close  to  1,  the  relative  error  and  the  relative  residual 


will  be  close  in  magnitude.  If  the  condition  number  is  large,  the  relative  error  could  be 
many  times  as  large  as  the  relative  residual.  The  condition  number  depends  on  the 
particular  vector  norm  that  is  used.  For  example,  for  the  2-norm  we  find  that  cond(A)  is 
simply  the  largest  singular  value  of  A  divided  by  the  smallest  singular  value  of  A  where 
matrix  A  is  not  necessarily  a  square  matrix.  The  proof  can  be  obtained  from  Leon  [42]. 
Assuming  next  that  8A  *  0  but  8b  =  0  in  Equation  4-14,  one  gets 

x  +  8x  =  (A  +  5A)~]  b  (4-24) 

Substitute  x  =  A~1b  into  Equation  4-24  and  rearrange  the  equation,  one  obtains 


8x  =  ((A  +  SAV]  -  A'1  )b 


(4-25) 


Use  the  matrix  identity 


B  1  -  A    =  A~l(A-B)B 


-i 


(4-26) 


One  gets 


x  =  A~l8A(A  +  8A)~]b 


(4-27) 


which  equals 


8x  =  A~*8A(x  +  8x) 


(4-28) 


Thus,  we  have 
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||<5x||<||a-J||m|||x  +  (5x| 


(4-29) 


or 


IN 


<  A  A 


-i 


INI 


(4-30) 


\x  +  8x 


\\4 


Similar  to  Equation  4-17  estimate  error  Sx  is  bounded  by  the  matrix  error  ratio  ||<5A||  /  ||A|| 
multiplied  by  the  condition  number.  Therefore,  the  condition  number  provides  the  upper 
boundary  of  the  estimation  error  in  relation  to  the  matrix  error  induced  by  the  parameter 
errors  and  measurement  noise. 

Most  the  least-squares  problems  are  overdetermined  problems  which  means  that 
the  matrix  is  not  a  square  matrix  but  a  rectangular  matrix.  The  condition  number  still  is 
applicable  for  this  case.  One  can  treat  the  matrix  inverse  of  A  as  a  Pseudoinverse  matrix 
which  is  A+  =  QL'UT  so  that  the  condition  number  derivation  process  from  Equation  4- 
14  to  4-30  is  still  valid. 


CHAPTER  5 
CALIBRATION  SIMULATION 


In  the  previous  chapter  we  introduced  three  important  calibration  steps  which  give 
us  the  basic  calibration  concept.  We  also  presented  the  observation  indices  that  provide 
us  a  reference  for  choosing  proper  HCMM  poses.  After  outlining  these  fundamentals, 
several  questions  then  are  raised:  How  best  to  select  measurement  configurations  or  what 
is  the  best  measurement  strategy?  How  many  poses  do  we  need  to  reach  the  desired 
accuracy?  What  is  the  effect  of  measurement  accuracy  and  noise?  Will  the  quality  of  the 
initial  estimate  of  parameters  be  a  factor  in  calibration  accuracy?  In  order  to  obtain  these 
answers  a  series  of  computer-based  mathematical  simulations  have  been  performed.  A 
simplified  HCMM  model  is  shown  in  Figure  5-1.  Ten  kinematic  parameters  which 
determine  the  initial  configuration  are  shown  in  Table  5-1.  The  maximum  and  minimum 
strut  length  are  52  and  32  inches,  respectively. 


Z  Axis 


(0,0.0) 


Figure  5-1  The  HCMM  model 
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Table  5-1  Initial  kinematic  parameters 

unit:  inch 


ai 

a2 

a3 

a4 

as 

a& 

43.10642755 

43.10642755 

43.10642755 

43.10642755 

43.10642755 

43.10642755 

r 

b 

h 

Lc 

68.233 

34.1165 

59.09151137 
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5.1  Simulation  Setup  and  Environment 


The  flowchart  for  the  calibration  simulations  is  shown  in  Figure  5-2.  The 
approach  taken  during  the  course  of  the  identification  simulations  was  as  follows.  The 
measurement  strategy  and  number  of  poses  were  decided  first.  The  initial  HCMM 
configuration  was  given  as  shown  in  Table  5-1.  The  center  rod  motion  then  was  created 
by  using  the  observation  strategy  (see  Section  5-2).  The  HCMM  reverse  kinematic  model 
then  was  applied  to  these  center  rod  positions  to  generate  the  strut  length  change  data, 
followed  by  a  check  of  the  strut  displacements  to  ensure  all  the  motions  are  within  the 
working  space  of  the  HCMM.  The  strut  length  changes  and  actual  center  rod  position 
coordinates  for  each  poses  then  were  saved  in  a  file,  with  14  place  decimal  accuracy.  A 
uniformly  distributed  random  noise  of  measurement  was  superimposed  on  the  strut  length 
changes  before  saving  to  the  data  file. 

The  programs  to  generate  measurement  data  and  to  implement  the  least-squares 
identification  method  were  written  and  ran  in  MATLAB  on  a  Windows  95  based  personal 
computer.  When  the  calibration  procedure  was  completed,  the  lower  sphere  position  at 
each  pose  was  estimated  by  applying  the  identified  kinematic  parameters  and  the  data  of 
the  stored  strut  length  changes  in  the  forward  kinematic  model.  The  lower  sphere 
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Observation  Strategy 


Number  of  poses 
strut  motion  range 


1 


Given  Center  rod  position 
and  orientation 

Pup,  Plow 


▼ 


compute  the  strut 
length  change  S 


nominal 
parameters 


1 

r 

f—  ► 

Least-Square  Parameter 
identification 

Measurement  noise 


T 

identified  parameters 
condition  number  of  Jacobian 


Figure  5-2  Flowchart  for  calibration  simulation 
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position  errors  were  determined  by  comparing  these  estimated  lower  sphere  positions  and 
the  stored  exact  lower  sphere  positions. 

The  symbols  used  in  the  simulation  results  table  are  as  follows.  The  number  of 
poses  is  "m."  The  error  value  used  for  the  ten  nominal  parameters  in  the  HCMM  model 
is  "8A."  The  measurement  uncertainty  is  represented  by  Un  which  is  assumed  as 
uniform  distribution  with  zero  mean.  The  measurement  strategy  types  are  "Strg."  After 
the  calibration  procedure  is  completed,  the  identified  parameter  errors  Sau     &6,  Sr,  5b, 
5h,  and  8Lc  are  represented  as  "&4rms"  which  is  found  by  taking  the  root  mean  square  of 
all  ten  parameter  errors.  The  charged  CPU  time  for  execution  is  "rcpu."  The  condition 
number  of  the  identification  Jacobian  is  "Cond(J)." 

In  our  iterative  least-squares  algorithm  each  entity  in  the  Jacobian  matrix  was 
calculated  analytically  by  taking  the  partial  derivative  of  objective  function  Equation  3- 
25.  The  iterative  procedure  was  terminated  when  the  maximum  value  in  the  parameter 
update  vector,  AA,  became  less  than  10"12  inch. 

5.2  Measurement  Strategy 

Several  measurement  strategies  have  been  examined  through  the  computer 
simulation.  The  number  of  poses  was  chosen  as  m=200  for  all  the  strategies.  A  noise 
level  Un=\Q\a  inch  was  superimposed  on  the  length  change  before  that  data  was  stored. 
A  initial  guess  error  8A=0.5  inch  was  added  to  each  of  the  exact  parameter  values. 
The  first  strategy  type  is  called  "straight  line  motion."  The  center  rod  was  moved  along  a 
straight  line  with  a  unit  vector  (0.6674,  0.6674,  0.3304)  and  was  kept  in  constant 
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orientation  during  the  movement.  The  strut  length  changes  were  stored  at  every  0.075 
inch  interval  along  the  center  rod  travel.  The  changes  in  length  are  shown  in  Figure  5-3. 

A  200x10  Jacobian  matrix  was  formed  by  the  pose  set .  Since  there  is  a  linear 
relationship  between  the  HCMM  configurations,  this  measurement  strategy  is  not 
sufficient  to  create  a  Jacobian  matrix  where  the  rank  of  J  is  less  than  10.  It  is  not 
identifiable  in  the  computer  simulation. 


10  30 

Figure  5-3  Straight  line  motion 

The  second  strategy  is  called  "2D  circular  motion,"  in  which  the  locus  of  the 
center  rod  forms  a  10  inch  diameter  cylinder.  Meanwhile  the  center  rod  was  kept  at  the 
same  z  height  and  parallel  to  the  z  axis.  The  strategy  is  shown  in  Figure  5-4.  The  same 
input  conditions  as  strategy  1  were  used  in  the  computer  simulation. 
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10  25 


Figure  5-4  2D  circular  motion 

Since  the  upper  and  lower  spheres  only  move  in  the  x  and  y  direction,  some 
kinematic  parameters  will  be  linearly  dependent  on  others.  Thus  the  rank  of  the  Jacobian 
matrix  from  this  measurement  strategy  would  be  non-full.  Consequently  the  parameters 
are  not  identifiable  by  using  this  strategy. 

The  first  two  strategies  only  involve  linear  or  planar  motion.  Even  the  200  poses 
used  could  not  provide  enough  information  for  identification.  Thus  a  squared  sine  wave 
motion  in  the  z  direction  was  added  to  the  2D  circular  motion,  where  the  center  rod 

moved  up  and  down  following  the  curve  z,  =  13  X  (sin(0,.  ))2 .  Therefore  this  third 
measurement  strategy  is  called  "3D  sine  wave  circular  motion",  shown  in  Figure  5-5. 
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Figure  5-5  3D  sine  wave  circular  motion 

The  same  input  conditions  as  strategy  1  were  used  in  this  computer  simulation. 
The  identified  parameter  error  dAms  is  0.00238  inch.  The  parameter  error  is  still  so  big 
that  this  strategy  would  not  be  adopted. 

In  the  first  three  strategies  the  center  rod  was  kept  parallel  to  the  z  axis.  The  forth 
measurement  strategy  is  designed  to  change  the  center  rod  orientation  during  data 
collection.  The  center  rod  was  tilted  45  degrees  from  z  axis.  Then  it  followed  the  same 
motion  type  as  the  third  strategy.  Therefore,  the  center  rod  orientation  will  vary  from 
position  to  position.  This  motion  is  then  called  "  tilted  3D  sine  wave  circular  motion,"  as 
shown  in  Figure  5-6.  The  same  input  conditions  were  used  here  as  in  the  other  three 
strategies.  The  identified  parameter  error  SA^  is  0.00210  inch,  which  is  still  in  the  same 
error  level  as  the  third  strategy. 
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Figure  5-6  Tilted  3D  sine  wave  circular  wave  motion 


The  last  measurement  strategy  was  called  "random  motion  strategy"  as  shown  in 
Figure  5-7.  Each  time  the  center  rod  was  moved  to  a  different  position  and  orientation 
randomly.  Since  there  is  no  clear  relationship  between  the  200  configurations,  this 
measurement  strategy  provided  sufficient  exciting  for  the  identification  matrix.  This  can 
be  seen  from  the  calibration  result  which  identified  the  parameters  with  error  as  of 
=  71|i  inch.  The  result  also  shows  a  smaller  condition  number  compared  with  the  other 
four  strategies.  The  comparison  of  all  the  five  strategies  is  listed  in  Table  5-2. 


Table  5-2  Comparison  of  measurement  strategies 


Strg  1 

Strg  2 

Strg  3 

Strg  4 

Strg  5 

SArms 

not  identifiable 

not  identifiable 

2378n  inch 

2100U  inch 

71(1  inch 

cond(J) 

ill  condition 

ill  condition 

1 1540 

12917 

430 

m=200  noise:- 1  On- 1  Onin        initial  8A=  0.5" 
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4(K 


0  20 

Figure  5-7  Random  motion 

5.3  The  Effect  of  Linear  Displacement  Measurement  Error  on  Sphere  Position 

As  we  know  the  linear  displacement  measurement  in  a  strut  is  performed  by  the 
laser  interferometer  system.  In  this  section  we  will  investigate  how  much  the  linear 
displacement  measurement  error  will  affect  on  the  accuracy  of  sphere  position  if  the  true 
initial  kinematic  parameter  values  are  known.  In  keeping  with  the  conventions  of 
mechanical  metrology,  the  measurement  error  should  be  treated  as  an  unknown  value.  It 
should  be  described  as  the  "measurement  uncertainty"  instead  of  using  the  word  "error." 
The  measurement  uncertainty  could  exhibit  any  distribution  pattern  including  normal 
distribution  or  uniform  distribution.  If  no  statistically  valid  description  could  be 
determined  for  this  measurement  uncertainty,  a  Type  B  uniform  measurement  error 
distribution  could  be  assumed  (see  Taylor  [42]).  Because  we  do  not  have  the  statistically 


56 

described  information  for  the  laser  interferometer  error,  the  Type  B  uniform  measurement 
error  is  assumed  and  will  be  used  in  our  remaining  simulation. 

In  the  simulation,  an  arbitrary  HCMM  configuration  was  chosen  from  within  the 
work  volume  of  the  machine.  If  the  linear  displacement  measurements  were  subjected  to 
uniformly  distributed  random  error  over  the  range  of  (-Un,  +Un),  the  absolute  strut  length 
error  will  be  the  same  as  the  linear  displacement  measurements  because  the  true  initial 
strut  length  was  assumed  to  be  known.  Thirty  thousand  different  configurations  were 
examined  ,then  the  sphere  position  error  in  each  configuration  were  calculated  and 
recorded.  The  averaged  position  error  in  each  of  the  X,Y,  and  Z  directions  was  calculated 
by  taking  the  mean  of  all  position  errors.  The  uncertainty  "ay"  associated  with  the 
position  error  is  found  by  taking  the  standard  deviation  of  the  position  errors.  The 
symbol  "i"  in  "ay"  represents  the  direction  (X.Y.or  Z),  and  "j"  represents  the  sphere 
(upper  or  lower).  Keeping  the  same  pose  but  doubling  the  uncertainty  level  the  position 
error  was  then  calculated  again.  Five  noise  levels  were  tested.  It  is  interesting  to  find  that 
the  distribution  of  the  position  error  in  30000  configurations  resulted  in  a  normal 
distribution  in  each  direction  for  each  sphere.  The  mean  value  of  the  position  error  was 
on  the  order  of  10"8  inch  in  all  the  directions  which  is  so  small  that  position  error  can  be 
treated  as  an  unbiased  normal  distribution.  The  position  uncertainties  of  the  spheres 
under  different  levels  of  linear  displacement  measurement  error  are  listed  and  plotted  in 
Table  5-3  and  Figure  5-8.  It  was  found  that  the  position  accuracy  is  linearly  dependent  on 
the  linear  displacement  measurement  uncertainty.  It  also  can  be  observed  that  positional 
uncertainty  in  the  X  or  Y  direction  is  approximately  half  of  the  measurement  uncertainty 
level,  and  in  the  Z  direction  has  the  same  order  of  magnitude  as  the  measurement 
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uncertainty  level.  The  small  residual  variances  in  the  position  error  when  Un  is  0  are 
attributed  to  computational  error  in  the  computer  simulation. 

Table  5-3  The  effect  of  strut  length  error  on  sphere  position 


Un 
(inch) 

<7xu 
(inch) 

(Tyil 
(inch) 

<3zu 
(inch) 

Cxi 
(inch) 

ayi 

(inch) 

azi 

(inch) 

0 

3.9xl015 

3.6xl015 

1.2xl014 

3.7xl015 

3.9xl0"15 

1.2xl014 

-5u~5u 

2.6u 

2.6u 

5.4u 

2.6ja 

2.6u 

5.6u 

-10u~  lOu 

5.2(1 

5.2n 

10.6)0. 

5.2u 

5.2u 

11. in 

-20(i  -  20(i 

10.2n 

10.2u 

22.0|i 

10.2u 

10.2u 

22.4u 

-40u  ~  40u 

20.8|i 

20.6(i 

44.6|i 

20.4u 

20.6(1 

45.3(1 

Position  uncertainty  (n  inch) 
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Figure  5-8  The  effect  of  strut  displacement  measurement  error  on  sphere  position 


5.4  The  Effect  of  Measurement  Uncertainty  on  Calibration  Results 


In  the  previous  section  we  assumed  that  the  true  initial  kinematic  parameter  values 
are  known,  and,  discussed  the  influence  of  the  strut  displacement  measurement 
uncertainty  on  the  positioning  uncertainty  of  the  spheres.  However,  in  reality  the  exact 
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values  of  the  kinematic  parameters  are  unknown  and,  therefore,  the  exact  sphere  position 
is  impossible  to  determine.  Thus  an  identification  procedure  is  needed  to  provided  the 
most  accurate  parameter  value.  The  following  test  is  designed  to  determine  the  effect  of 
the  strut  displacement  measurement  noise  on  the  kinematic  parameter  identification 
process.  To  control  the  experimental  factor,  the  poses  used  in  parameter  calibration 
procedure  will  be  fixed.  The  only  variable  left  is  the  magnitude  of  the  measurement 
uncertainty  level  Un.  Thirty  poses  were  generated  by  using  strategy  5.  Uniformly 
distributed  random  noise  with  magnitude  Un  was  then  superimposed  on  the  strut  length 
changes  to  simulate  the  measurement  uncertainty.  The  approximate  or  nominal 
parameter  values  may  be  obtained  by  simply  using  a  tape  measure.  In  this  simulation  a 
0.5  inch  initial  parameter  error  "5A"  was  added  to  each  of  the  ten  true  parameters  of  the 
HCMM  as  the  nominal  parameter  values  for  least-squares  algorithm.  After  the 
calibration  procedure  was  completed,  the  root  mean  square  value  of  the  ten  kinematic 
parameter  errors  was  stored  as  "SAi^"  The  position  errors  in  upper  and  lower  spheres 
also  were  recorded  "8P." 

In  order  to  evaluate  the  effect  of  noise  magnitude  on  sphere  position  estimation, 
five  hundred  simulations  were  performed.  Each  simulation  run  used  the  same  poses  and 
noise  magnitude.  The  results  were  different  for  each  simulation  because  the  randomly 
generated  noise  was  different  each  time.  The  simulation  flowchart  is  shown  in  Figure  5-9. 

The  relationship  between  strut  displacement  measurement  uncertainty  and  the 
calibration  results  then  could  be  addressed  by  statistical  means.  From  Table  5-4  and 
Figure  5-10  it  can  be  observed  that  when  the  measurement  uncertainty  level  was  doubled, 
the  calibration  errors  also  doubled.  This  test  also  shows  that  the  parameter  values  could 
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Given  30  HCMM  poses  by 
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Figure  5-9  Flowchart  for  finding  the  relation  between  the  strut  displacement 
measurement  uncertainty  and  calibration  results 
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Table  5-4  Calibration  results  subjected  to  different  strut  displacement  measurement 
uncertainty  level 

Unit:  inch 


meas.  Uncert. 
Un 

mean  of 

mean  of 
cond(J) 

Ocond(J) 

0 

2.31x10° 

1.85xlO"20 

453.53 

4.10xl0'3 

-5u~5(i 

74.7U 

31.8(1 

453.53 

4.98xl0-3 

-10n~10u 

152(1 

66.4(1 

453.53 

5.00x1 0"3 

-20(i~20|i 

302(1 

131U 

453.53 

7.63X10"3 

-40|i~40(i 

619U 

280fl 

453.53 

1.40X10"2 

Un 

azu 

0 

2.78xl0"13 

6.76X1014 

3.16xl0"14 

2.77X10"13 

6.94xl0"13 

3.08X1013 

-5u~5(i 

76u 

20.8U 

lOu 

77(1 

24.1(1 

9.5(1 

-10|i~10|i 

1 50(i 

42.8u 

19u 

154(1 

51.5(1 

18.2(1 

-20|i~20|i 

316u 

89.3(1 

40.7(1 

317(1 

100(1 

37.3(1 

-40(i~40(i 

638u 

178u 

82(1 

648(1 

200(1 

79.2(1 

unique  pose  set  m=30  N=500  times      initial  8A=  0.5" 


Figure  5-10  The  effect  of  the  strut  displacement  measurement  error  on  the 
calibration  results  (Unique  pose  set,  Cond=453,  5A=0.5",  N=500  simulations) 
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be  identified  to  nearly  their  true  values  if  no  measurement  error  was  applied.  But  once 
the  measurement  error  was  introduced,  it  became  impossible  to  obtained  the  exact 
parameter  values. 

By  using  the  same  calibration  poses  only  a  slight  increase  was  observed  in 
condition  number  when  the  measurement  uncertainty  went  up.  This  shows  that  the 
condition  number  was  mainly  determined  by  the  selected  poses.  It  also  can  be  observed 
from  Figure  5-10  that  both  of  the  upper  and  lower  spheres  in  the  Z-direction  have  the 
smallest  positional  uncertainty.  The  X-direction  has  the  largest  positional  uncertainty. 
Why  does  this  occurs?  To  explain  this  phenomena,  we  trace  back  to  the  calibration 
method.  The  objective  function  of  the  self-calibration  algorithm  is  to  find  a  set  of 
parameters  such  that  the  center  rod  length  has  the  minimum  deviation  in  all  the 
configurations.  After  we  performed  the  calibration,  the  least  squares  method  found  an 
optimal  set  of  parameters  that  made  the  deviation  in  the  center  rod  length  the  smallest. 
As  we  know,  the  center  rod  length  is  derived  from  a  set  of  equations,  Equation  3-2  to 
Equation  3-8.  All  the  ten  kinematic  parameters  are  necessary  to  compute  the  center  rod 
length  deviation.  However,  the  X-coordinate  is  the  first  equation  in  the  series  of 
equations  used  to  obtain  the  center  rod  length.  It  only  uses  four  of  the  ten  kinematic 
parameters  as  shown  in  Equation  3-2  and  Equation  3-5.  The  search  procedure  does  not 
find  the  four  parameters  which  yield  the  best  estimate  for  the  X  coordinate,  but  rather  the 
center  rod  length  with  the  least  deviation.  In  contrast,  the  Z  coordinate  is  the  final 
equation  in  the  series.  Six  parameters  are  needed  to  compute  each  Z-coordinate. 
Furthermore,  since  the  center  rod  orientation  never  deviates  more  than  30°  from  the 
vertical,  the  center  rod  length  is  more  strongly  dependent  on  the  Z  coordinate  of  the 
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spheres  than  it  is  on  the  X  and  Y  coordinates.  Therefore,  it  can  be  expected  that  Z  has  the 
smallest  uncertainty  level  after  calibration. 

It  should  be  noted  that  a  calibration  using  different  calibration  pose  sets  could 
produce  a  different  uncertainty  level  in  the  resulting  estimation  of  the  ten  kinematic 
parameters. 

5.5  The  Effect  of  Pose  Selections  on  Calibration  Results 

Each  pose  set  has  a  unique  sensitivity  level  for  the  measurement  noise,  such  that 
using  a  different  pose  set  would  yield  different  calibration  results.  To  determine  the  noise 
sensitivity  of  a  pose  set,  a  widely  used  index  is  the  condition  number  of  the  Jacobian 
matrix  J.  In  the  previous  simulation  the  identical  pose  set  was  used  to  find  the 
relationship  between  measurement  uncertainty  and  calibration  results.  In  order  to 
understand  the  effect  of  pose  set  selection  on  the  calibration  results,  a  thousand 
calibration  simulations  were  performed.  For  each  simulation  a  pose  set  consisting  of  30 
poses  was  generated  by  using  a  random  selection  strategy.  The  flowchart  of  this 
simulation  is  shown  in  Figure  5-11. 

The  calibration  results  for  Un=±10|i  inch  are  shown  in  Figure  5-12.  Each  circle 
point  in  Figure  5-12  represents  a  simulation  result  of  §Ams.  Each  pose  set  is  associated 
with  its  condition  number  which  represents  its  sensitivity  level  to  noise. 

One  can  observe  from  the  simulation  results  that  a  pose  set  with  smaller  condition 
number  doses  not  guarantee  that  the  calibration  results  will  be  better  than  the  pose  set 
with  higher  condition  number.  However,  comparing  the  upper  boundaries  of  condition 
numbers,  one  can  see  the  trend  that  the  pose  sets  with  smaller  condition  numbers  have 
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Given  30  HCMM  poses  by 
using  Strategy  5 
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Figure  5-1 1  Flowchart  of  calibration  simulation  designed  for 
understanding  the  effect  of  pose  selection  on  the  calibration  results 
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1.00E-03 


Condition  Number 

Figure  5-12  Calibration  results  and  the  coresponding  condition  number  for  each  pose  set 
(m=30  poses,  Un=±10^i  inch,  N=1000  simulations) 

smaller  upper  bound  in  position  errors.  If  we  draw  a  straight  line  above  all  the  points  in 
Figure  5-12,  we  can  see  the  trend.  A  more  detailed  discussion  of  this  issue  will  described 
in  Section  5-7. 

From  the  above  simulation  we  found  that  SA^  varies  from  10.4  |i  inch 
(cond=53 1)  to  916|X  inch  (cond=1430).  However,  for  a  randomly  selected  pose  set  the 
calibration  result  has  a  large  possibility  of  being  located  in  the  region  of  the  mean  value 
from  1000  simulation  results.  Thus  we  choose  to  use  the  mean  value  and  its  standard 
deviation  of  the  1000  simulation  results  to  represent  the  overall  calibration  results.  This 
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simulation  was  repeated  for  different  levels  of  measurement  uncertainty.  The  results  are 
listed  in  Table  5-5. 


Table  5-5  The  relation  between  the  strut  displacement  measurement  uncertainty 
and  calibration  results 

Unit:  inch 


meas.  uncert. 
Un 

mean  of 
SAnm 

mean  of 
Cond(J) 

^Cond(J) 

0 

3.24xl0"13 

2.37X10-20 

713 

150.4 

-5(1  -  5u 

104(1 

52.9m- 

712 

154.4 

-10(1  ~10u 

206|i 

101(1 

710 

148.8 

-20(1  ~  20(i 

426u 

206(1 

715 

146.2 

-40u  ~  40u 

817(1 

395(1 

707 

148.3 

m=30  poses 

N=  1000  times 

initial.&4  =  0.5 

This  result  is  similar  to  the  previous  section's  result  in  that  the  overall  calibration 
results  are  linearly  dependent  on  the  measurement  uncertainty.  These  results  also  show 
that  in  the  absence  of  measurement  uncertainty,  the  parameters  could  be  exactly  identified 
as  long  as  its  Jacobian  matrix  is  not  ill  conditioned. 

5.6  The  Effect  of  Number  of  Poses  on  the  Calibration  Results 

To  investigate  the  importance  of  the  number  of  poses  to  the  calibration  results, 
these  calibrations  are  separated  into  five  groups  with  m=30,  50,  100,  200,  and  300  poses 
respectively.  The  measurement  uncertainty  is  assumed  to  be  ±10  |iin  for  each  of  the  five 
cases,  and  the  initial  nominal  parameter  error  is  0.5  inch  for  all  parameters.  For  the  cases 
of  m=30,  50  and  100  poses  a  total  of  2000  simulations  for  each  group  were  performed. 
The  poses  in  each  simulation  are  regenerated  by  using  strategy  5.  For  the  cases  of  m=200 
and  300,  a  total  of  1000  simulations  were  performed. 

The  results  shown  in  Table  5-6  demonstrate  that  with  an  increasing  number  of 
poses  in  the  calibration  procedure,  a  decreasing  condition  number  will  result.  Thus  better 
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calibration  results  can  be  expected.  It  also  can  be  observed  that  when  the  number  of 
poses  is  more  than  200,  the  improvement  in  calibration  result  is  minimal. 

Table  5-6  The  Effect  of  number  of  poses 


unit:  inch 


no.  of 
poses  (m) 

mean  of 

^SArms 

mean  of 
cond(J) 

^  Cond(J) 

CPU 
time(sec) 

30 

206u 

101m 

710.1 

148.8 

0.86 

50 

143p. 

66.5jll 

566.8 

87.6 

1.56 

100 

100m 

44.0m 

481.3 

58.1 

2.90 

200 

71.6m 

32.2m 

436.4 

41 

5.86 

300 

59.2m 

26.1m 

417.0 

35.5 

9.31 

Strategy:  5         Un=±10nin  initial  8A=0.5"    CPU:  Pentium  200 


The  diamond  line  in  Figure  5-13  shows  the  mean  of  SArmS  from  the  calibration 
simulation  results,  and  the  square  line  shows  the  error  amount  after  adding  the  standard 
deviation  to  the  mean  position  error.  If  the  measurement  uncertainty  is  within  ±10  Min 
and  the  pose  set  is  randomly  selected,  it  can  be  expected  that  the  calibration  result  in 
SArms  will  have  a  68.27%  chance  of  falling  below  the  square  line.  The  maximum  8Arms 
from  each  of  the  2000  simulations  were  recorded  as  the  local  upper  bound.  This  bound 
gives  us  a  general  idea  that  most  of  the  calibration  results  will  not  exceed  these  values. 

Figure  5-14  shows  the  relation  between  the  condition  number  and  the  number  of 
poses  used  in  the  calibration.  The  more  poses  which  are  used  in  the  calibration,  the 
smaller  the  deviation  in  the  condition  number.  Thus,  the  condition  number  is  gradually 
losing  its  ability  in  determine  the  goodness  of  a  pose  set  when  the  number  of  poses  is 
greater  than  200.  This  phenomena  will  be  discussed  in  the  next  section. 
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Figure  5-13  Parameter  error  vs  number  of  poses  (Un=±10m  inch) 
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Figure  5-14  Condition  Number  vs.  Number  of  Poses 
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5.7  The  Sensitivity  Analysis  on  Condition  Number 

In  previous  sections  we  observed  from  the  simulation  results  that  a  pose  set  with 
smaller  condition  number  does  not  guarantee  that  the  calibration  result  will  be  better  than 
the  pose  set  with  higher  condition  number.  We  also  observed  that  the  more  poses  which 
were  used  by  the  calibration  routine,  the  smaller  the  deviation  in  the  condition  number. 
In  this  section  we  will  investigate  how  the  condition  number  can  help  us  in  selecting  a 
pose  set  that  will  have  a  better  chance  to  obtain  good  calibration  results,  and  the  limit  to 
the  condition  number's  usefulness. 

According  to  the  magnitude  of  the  condition  number,  we  selected  several  groups 
of  pose  sets.  Each  group  contains  two  to  three  pose  sets,  and  each  pose  set  in  the  same 
group  contains  30  poses  with  condition  number  which  varies  in  range  ±5.  For  each  pose 
set  300  calibration  simulations  were  performed  as  shown  in  Figure  5-9.  In  each 
simulation,  a  uniformly  distributed  measurement  uncertainty  of  ±10  pin  was 
superimposed  on  the  strut  length  changes,  and  the  initial  parameter  error  was  0.5  inch  for 
all  the  ten  kinematic  parameters.  Each  of  the  calibration  simulation  results  was 
represented  as  a  circle  as  shown  in  Figure  5-15.  The  calibration  results  in  each  group 
were  represented  in  a  statistical  manner  as  listed  in  Table  5-7  and  plotted  in  Figure  5-16. 

To  verify  the  accuracy  of  the  statistical  representation  for  condition  number,  we 
performed  ten  thousand  calibration  simulations.  In  each  simulation,  we  arbitrarily 
selected  30  poses  and  applied  the  same  noise  level  and  initial  parameter  errors  as  in  the 
above  simulations.  The  accuracy  of  the  statistical  representation  was  verified  by  plotting 
the  ten  thousand  simulation  results  on  top  of  Figure  5-16  as  shown  in  Figure  5-17. 
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Figure  5-15  Test  of  condition  number 


Table  5-7  Statistical  representation  for  condition  number 


unit:  |J.  inch 


condition 

mean  of 

std.  dev. 

minimum 

maximum 

Number 

SArms 

SArms 

388 

145 

59 

29.0 

281 

504 

160 

70 

24.2 

391 

602 

169 

85 

18.0 

461 

699 

192 

96 

21.4 

537 

805 

220 

103 

32.5 

728 

1003 

262 

144 

35.0 

859 

1098 

278 

142 

16.1 

894 

1203 

296 

165 

20.5 

1060 

m=30  initial  8A=0.5" 
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Figure  5-16  The  statistical  representation  for  condition  number 
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Figure  5-17  Verification  of  the  statistical  representation  for  condition  number 
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Figure  5-17  to  Figure  5-19  show  the  statistical  representations  for  the  calibration 
simulation  results  by  using  different  numbers  of  poses  in  the  calibration.  It  can  be 
observed  that  more  poses  used  in  the  calibration  resulting  in  a  smaller  scattering  zone  for 
the  calibration  results.  The  comparison  of  the  condition  number  sensitivity  is  shown  in 
Figure  5-20.  Each  line  in  Figure  5-20  is  the  mean  statistical  line  for  m=30,50,100,  and 
200.  At  m=30,  a  pose  set  with  a  smaller  condition  number  will  statistically  have  better 
calibration  results  than  a  pose  set  with  a  higher  condition  number.  But,  for  the  case  of 
m=200,  a  pose  set  with  a  smaller  condition  number  does  not  mean  it  will  have  better 
calibration  results  than  a  pose  set  with  a  higher  condition  number.  Therefore,  we  realize 
that  condition  number  loses  its  sensitivity  when  more  poses  are  used  in  calibration. 

To  choose  an  optimal  pose  set  with  30  poses,  we  performed  ten  thousand 
calibration  simulations  and  used  the  condition  number  as  our  searching  index.  The  best 
pose  set  we  found  has  a  condition  number  of  388.  Another  best  pose  set  with  50  poses 
also  was  found  using  the  same  method  which  has  a  condition  number  of  347.  These  two 
best  pose  sets  are  listed  in  the  Appendix  B. 

5.8  The  Effect  of  Nominal  Value  on  the  Calibration  Results 

To  facilitate  a  numerical  solution  method,  the  nonlinear  least-squares  equation  is 
linearized  at  its  nominal  values.  Therefore,  it  is  necessary  to  understand  the  effect  of  the 
nominal  parameter  error  on  the  calibration  results.  Three  pose  sets  were  tested.  Each 
pose  set  includes  thirty  poses  selected  by  strategy  5.  The  nominal  parameter  errors  were 
increased  from  zero  to  more  than  three  inches.  Table  5-8  shows  the  calibration  results 
were  identical  as  long  as  the  least-squares  objective  function  does  not  diverge.  The 
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Figure  5-17  Statistical  representation  for  calibration  simulation  results  (m=50) 
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Figure  5-18  Statistical  representation  for  calibration  simulation  results  (m=100) 
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Figure  5-19  Statistical  representation  for  calibration  simulation  results  (m=200) 
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Figure  5-20  Sensitivity  of  the  condition  number 
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divergence  will  happen  when  the  parameter  errors  are  large  enough  to  cause  Equation  3-4 
or  Equation  3-7  to  become  a  complex  number.  The  divergence  can  be  avoided  by 
choosing  a  more  accurate  nominal  parameter  value  or  by  selecting  a  different  pose  in 
which  the  complex  number  does  not  appear  in  the  objective  function.  However,  in  all 
cases  where  the  objective  function  does  in  fact  converge,  the  calibration  result  will  not  be 
affected  by  the  initial  nominal  parameter  error. 

Table  5-8  The  effect  of  initial  nominal  value  error 


unit:  in 


initial  nominal 

pose  set#l 

pose  set  #2 

pose  set  #3 

parameter  error 

SArms 

SAjms 

0 

173.8)1 

144.2m 

193.6m 

1 

173.8n 

144.2m 

193.6m 

1.5 

173.8m. 

144.2m 

193.6m 

2 

173.8m 

144.2m 

193.6m 

2.5 

173.8m 

144.2m 

193.6m 

3 

173.8m 

diverge 

diverge 

above  3 

diverge 

diverge 

diverge 

m=30  Un=+10|lin  strategy:  5 


5.9  The  Comparison  of  Identification  Methods 

As  mentioned  earlier,  many  methods  can  be  used  to  solve  for  the  nonlinear  least- 
squares  problem.  The  two  most  widely  used  methods  are  the  linearized  least-squares 
method  and  the  Levenberg-Marquardt  method  (LM  method)  which  is  a  non-linear 
method,  or  a  modification  of  Newton's  method.  The  simulations  were  run  in  MATLAB 
under  Windows  95  on  a  Pentium  200  personal  computer.  Table  5-9  shows  the 
comparison  of  these  two  solution  methods.  Both  the  linearized  and  the  non-linear 
method  (LM  method)  gave  almost  identical  calibration  results  under  the  same  conditions. 
This  table  takes  twelve  cases  of  increasing  numbers  of  poses  and  compares  the  CPU  time 


taken  by  the  two  methods.  The  linearized  method  is  seven  to  twelve  times  faster  than  the 
non-linear  LM  method.  At  200  poses  the  linear  method  takes  about  six  seconds,  while 
the  LM  method  takes  more  than  48  seconds  of  CPU  time. 


Table  5-9  The  comparison  of  two  identification  methods 


Linearized  method 

LM  method 

^cpu  (sec) 

SArms  (Ltin) 

t^py  (sec) 

SArms  (llin) 

test  1 

m=30        test  2 
test  3 

0.86 

193.606 

7.8 

193.606 

0.82 

188.715 

7.8 

188.727 

0.88 

253.861 

8.41 

253.861 

test  1 

m=50         test  2 
test  3 

1.76 

101.831 

13.35 

101.843 

1.15 

88.579 

10.82 

88.582 

1.16 

107.816 

14.01 

107.822 

test  1 

m=100       test  2 
test  3 

2.69 

34.585 

26.47 

34.585 

2.67 

83.623 

19.61 

83.619 

3.06 

47.473 

21.26 

98.731 

test  1 

m=200        test  2 
test  3 

5.43 

39.840 

55.3 

39.840 

5.64 

33.7 

54.82 

33.698 

6.21 

44.089 

48.17 

44.087 

5.10  Discussion 

The  kinematic  parameters  can  be  exactly  identified  only  in  the  theoretical  case 
which  has  no  displacement  measurement  error  and  modeling  error.  When  measurement 
uncertainty  exist  it  is  impossible  for  one  to  determine  the  true  parameter  values. 

A  good  measurement  strategy  will  improve  the  quality  of  the  parametric 
identification  results.  In  some  cases,  the  parameters  are  unidentifiable  from  a  given  pose 
set,  such  as  the  straight  line  or  circular  motion.  The  best  measurement  strategy  of  those 
tested  is  random  selection.  Furthermore,  it  has  been  found  that  not  all  of  the  pose  sets 
selected  from  the  random  strategy  obtained  acceptable  calibration  results.  A  good 
calibration  result  should  be  obtained  when  the  pose  set  is  least  sensitive  to  strut 
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displacement  measurement  noise.  From  simulation,  it  shows  that  the  calibration  results 
are  linearly  decreased  with  the  increase  of  the  strut  displacement  measurement 
uncertainty. 

It  has  been  shown  that  increasing  the  number  of  poses  and  reducing  the 
measurement  error  are  two  methods  to  improve  the  calibration  quality.  In  general,  the 
quality  of  initial  nominal  parameter  guesses  will  not  affect  the  identification  result.  When 
the  number  of  poses  are  limited  to  less  than  200,  a  significant  factor  that  will  improve  the 
calibration  quality  is  the  pose  set  selection.  Thus,  the  condition  number  of  the 
identification  matrix  can  be  used  as  an  important  indicator  during  the  selection  of  the 
poses.  The  condition  number  computed  from  the  experiment  can  be  used  to  refer  back  to 
the  simulation  tables  and  estimate  the  quality  of  the  pose  set,  and  suggest  a  level  of 
confidence  in  parameter  identification. 


CHAPTER  6 
THERMAL  EFFECTS 

From  the  previous  chapter,  it  is  known  that  the  self-calibration  method  is  very 
sensitive  to  the  measurement  error  of  the  strut  linear  displacement.  When  the 
environmental  temperature  varies,  the  components  in  the  HCMM  will  exhibit  a  change  in 
their  dimensions.  Though  the  whole  calibration  procedure  is  assumed  to  be  completed  in 
a  very  short  time  period,  the  thermal  effect  inevitably  will  have  some  affect  on  the 
accuracy  of  the  linear  displacement  measurement  in  each  strut.  Thus,  to  investigate  the 
thermally  induced  error  of  the  linear  displacement  measurement  in  a  single  strut  is  one 
objective  of  this  chapter.  The  self-calibration  method  is  based  on  the  center  rod 
remaining  a  constant  length  during  the  calibration  operation  period.  Since  the  center  rod 
will  experience  thermally  induced  changes,  qualification  of  the  thermal  effect  on  the 
center  rod  is  a  second  objective  of  this  chapter. 

The  coordinates  of  a  probe  are  obtained  by  knowing  the  absolute  lengths  of  the  six 
struts  and  the  three  base  sphere  distances.  If  these  lengths  are  not  corrected  when  the 
thermal  effect  occurs,  the  length  errors  will  propagate  through  the  probe  coordinate 
calculations  and  cause  errors  in  the  measured  position  of  the  probe.  Thus,  understanding 
the  thermal  deformation  in  each  component  is  another  objective  of  this  chapter. 
Furthermore,  thermal  effects  may  increase  over  time.  Thus,  to  determine  the  time  period 
beyond  which  the  HCMM  needs  to  be  re-calibrated,  when  the  thermal  uncertainty  may 
have  exceeded  the  allowable  tolerance  is  included. 
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Since  the  actuator  is  used  to  control  the  motion  of  the  strut,  it  must  generate  some 
heat  during  its  operation.  To  investigate  if  the  generated  heat  in  a  motor  will  cause  a 
significant  effect  in  the  measurement  error  of  the  strut  length  is  the  last  subject  in  this 
chapter. 

Assumptions  must  be  made  before  any  thermal  analysis  is  undertaken.  As  we 
know  the  net  thermal  effect  can  be  represented  as  the  combination  of  effects  from  three 
distinct  mode  of  heat  transfer:  conduction,  convection,  and  radiation.  We  do  not  expect 
there  to  be  any  high  temperature  heat  sources  nearby  or  in  contact  with  the  HCMM. 
Therefore,  for  simplicity  in  the  thermal  analysis,  we  assume  the  effect  of  radiation 
between  the  environment  and  the  HCMM  can  be  ignored.  The  only  thermal  effect  on  the 
HCMM  considered  in  this  chapter  is  the  natural  convection  between  the  ambient 
atmosphere  and  the  surface  of  the  HCMM.  We  also  assume  that  the  ambient  temperature 
is  uniformly  distributed  throughout  the  environment,  such  that  the  ambient  temperature  is 
considered  to  be  a  function  of  time  only. 

6. 1  Thermally  Sensitive  Components 

The  absolute  length  of  a  strut  is  measured  from  the  base  sphere  center  to  the 
moving  sphere  center.  Strut  components  include  two  reference  spheres,  shoulder, 
actuator,  ball  screw,  and  the  pin  followers.  The  components  and  the  associated 
dimensions  of  a  strut  are  shown  in  Figure  6-1.  In  a  strut,  the  laser  interferometer  system 
will  detect  a  length  change  from  its  emitter  to  the  retroreflector.  However,  the  thermal 
growth  in  the  components  which  is  not  covered  by  the  laser  path  will  not  be  detected. 
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This  unmeasured  thermal  growth,  then,  will  be  treated  as  a  noise  or  measurement  error  in 
the  strut  length. 


Base  Spnere 

Pin  follower  0-25"D 


Shoulder  0D=2.5" 
 10=2.22" 


Ball  Screw  OD-r      /       Upp.r  (Lower) 
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Retroreflector/ 
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sensttve 
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Loser  Path 

Thermal  growth  can  be  detected  by  laser  interferometer 


Absolute  strut  length 


thermal 
senstrve 
area 


Figure  6- 1  The  thermal  sensitive  components  in  a  strut 

The  center  rod,  as  shown  in  Figure  6-2  contains  three  different  material  elements: 
reference  spheres,  joint  holder,  and  fiberglass  tube.  When  the  ambient  temperature 
changes,  the  center  rod  will  change  its  length  so  that  the  constant  center  rod  length 
assumption  will  be  violated.  Understanding  the  thermal  response  of  the  current  center 
rod  design  will  help  us  in  designing  a  new  center  rod  if  the  redesign  is  necessary. 


=0 


.Ref.  Sphere  2"D 


Fiberglass  1"Dx29"L 


Aluminum 
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Figure  6-2  The  center  rod 


Because  the  base  spheres  are  fixed  on  the  triangle  frame  as  shown  in  Figure  6-3, 
any  thermal  growth  in  the  triangle  frame,  will  cause  the  distance  between  the  base  spheres 
to  change.  If  the  thermal  growth  is  not  corrected  from  the  nominal  values,  it  will 
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propagate  to  increase  the  uncertainty  of  the  calculated  probe  coordinates.  Symmetry 
prevents  the  triangle  frame  from  being  distorted  when  it  is  subjected  to  an  environmental 
temperature  change.  Thus,  equal  thermal  growth  rate  is  assumed  for  the  thermal  effect 
analysis. 


60 


3.75W  x  3.75"H 


60 


Base  sphere  length  68.33" 

Figure  6-3  The  triangle  frame 


The  supporting  frame  is  another  component  that  needs  to  be  considered.  From 
Figure  6-4,  one  can  observe  that  any  thermal  growth  in  the  columns  will  make  the  triangle 
frame  move  up  or  down.  Therefore  any  relative  position  changes  between  the  world 
coordinate  system  and  the  granite  table  will  be  directly  added  onto  the  measurement  error 
in  the  Z-direction,  but  not  through  the  error  propagation  procedure. 
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6.2  The  Thermal  Response  in  The  Sensitive  Components 

For  a  shop  floor  with  conventional  heating  or  cooling  systems,  the  ambient 
temperature  could  be  idealized  as  shown  in 

Figure  6-5.  The  temperature  variation  range  and  its  cycle  frequency  will  depend 
on  the  setup  of  the  heating  or  cooling  systems.  In  this  case,  it  is  important  to  know  the 
thermal  response  of  the  HCMM  components  under  typical  ambient  variation. 


Figure  6-5  The  temperature  variation  in  a  shop  floor 


Review  of  Heat  Transfer 

To  analyze  the  thermal  response,  according  to  Incropera  and  DeWitt  [43],  first  we 
check  the  Biot  number  of  a  solid  component  in  its  thermal  model.  The  Biot  number  is  a 
dimensionless  group  and  is  defined  as  shown  in  Equation  6.1. 


(6-1) 
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The  hc  is  the  thermal  convection  coefficient  between  the  air  and  the  surface  of  the 
solid  component.  The  Ls  is  the  characteristic  length  which  is  defined  as  the  ratio  of  the 
solid  volume  to  surface  area,  Lc=  V/As.  The  k  is  the  thermal  conductivity  of  the  solid. 
For  a  Biot  number  which  is  much  less  than  1.0,  the  thermal  model  can  be  treated  as  a 
lumped  system.  In  the  lumped  system  model  it  is  assumed  that  the  temperature 
distribution  within  the  solid  at  any  instant  is  sufficiently  uniform  that  the  temperature  of 
the  solid  can  be  considered  to  be  a  function  of  time  only.  For  a  the  Biot  number  less  0.1, 
the  assumption  of  uniform  temperature  distribution  has  an  error  less  than  5  percent. 

For  a  lumped  system  the  energy  equation  for  heat  transfer  in  the  solid  may  be 
stated  as  :  (rate  of  heat  flow  into  the  solid  of  volume  V  through  boundary  surfaces  A)  = 
(rate  of  increase  of  internal  energy  of  the  solid  of  volume  V).  By  writing  the  appropriate 
mathematical  expression  for  each  of  these  terms,  it  becomes 

hcAs[T„-T(t)]  =  pVc^p-  (6-2) 

where  p=  kg/m3  is  the  density  of  the  solid,  c=(J/(kg.  °C)  is  the  specific  heat  of  materials. 
V  is  its  volume,  As  is  its  surface  area.  Too  is  the  ambient  temperature. 

To  solve  Equation  6-2  we  must  specify  the  temperature  of  the  body  at  a  particular 
time.  If  we  assume  that  the  temperature  of  the  body  at  an  initial  time,  t=0,  is  known  to  be 
T0,  we  say  that  the  initial  condition  for  Equation  6-2  is  0O  =  T0-Too  at  t  =  0,  and  0(t)  = 
T(t)  -  T„„.  The  solution  to  Equation  6-2  is 

_  p-(hcAslpVc)t 

a      ~e  (0-3) 

O 

Equation  6-3  can  be  rewritten  as 
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M=e-«  (6.4) 

e. 

where  T  represents  thermal  time  constant  of  the  solid,  and  its  value  is 
pVc 

t=mT  (6"5) 

When  the  environmental  temperature  is  known,  the  temperature  in  the  solid 
component  can  be  estimated  by  using  Equation  6-4.  For  instance,  a  steel  ball  with  a 
uniform  temperature  of  T0  is  suddenly  moved  into  a  room  where  the  air  temperature  is 
well  maintained  at  Too  all  the  time.  The  time  required  for  the  steel  ball  to  rise  63.2%  of 
the  temperature  difference  between  the  air  and  steel  ball  initial  temperature  is  called  the 
time  constant,  as  shown  in  Figure  6-4(a).  The  temperature  of  the  steel  ball  at  time  t  can 
be  estimated  using  Equation  6-4,  and  can  be  written  as 

T(t)  =  T„  +  (T0-T„  )e~'H  (6-6) 


87 

For  an  environment  with  a  sinusoidally  varying  temperature,  the  temperature 
curve  can  be  discretized  by  dividing  the  time  into  several  small  steps.  During  each  time 
step  the  air  temperature  is  assumed  to  remain  constant,  as  shown  in  Figure  6-4(b).  Thus, 
Equation  6-6  can  be  applied  in  this  case  and  rewritten  as 

At 

Tn+l  =  Tair(n+l)  +  (Tn  ~  Tair(n^)  )e    %  (6-7) 

where  n  is  an  integer  starting  from  0.  At  is  the  time  interval  of  a  step. 

Once  the  temperature  of  the  solid  is  obtained,  the  thermal  growth  of  the  solid  can 
be  calculated  by  the  following  formula 

AL  =  ax  LxAT  (6-8) 
where  a  is  the  thermal  expansion  coefficient  of  the  solid,  L  is  the  length  of  the  solid  in 
the  thermal  growth  direction,  and  AT  is  the  temperature  change  of  the  solid  from  the 
initial  temperature. 

Lumped  System  and  Thermal  Time  Constant 

To  verify  if  the  lumped  system  model  is  applicable  on  the  HCMM's  thermally 
sensitive  components,  and  furthermore  to  obtain  the  time  constant  of  each  component, 
the  physical  properties  of  each  component  in  HCMM  have  to  be  obtained  first.  The 
listed  physical  property  of  a  material  in  different  reference  handbooks  show  different 
values  if  its  exact  chemical  formulation  is  not  available.  The  physical  properties  shown 
in  Table  6-1  are  the  average  values  taking  from  different  reference  resources  [43-46].  For 
instance,  for  stainless  steel,  the  uncertainty  range  is  about  2%  for  density  (p),  10%  for 
thermal  conductivity  (k),  8%  for  specific  heat  (c),  and  10%  for  the  thermal  expansion 
coefficient  a).  Here  we  take  the  average  from  the  available  data  in  different  resources. 
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The  uncertainty  level  for  aluminum  is  about  2%  for  density  (p),  20%  for  thermal 

conductivity  (k),  10%  for  specific  heat  (c),  and  2%  for  the  thermal  expansion  coefficient 

(a).  The  material  properties  for  fiberglass  tube  is  not  available  from  the  handbooks 

mentioned  in  above.  Its  thermal  expansion  coefficient  is  obtained  through  a  simple 

experiment.  A  small  section  of  fiberglass  tube  of  0.12616m  was  measured  at  25°C.  Then 

it  was  placed  in  a  freezer  for  a  long  time  such  that  its  temperature  was  close  to  the 

freezer's  temperature,  0°C.  The  fiberglass  tube  was  measured  again,  and  its  length  was 

0.12604m.  If  the  fiberglass  tube  is  assumed  to  have  the  same  temperature  as  the  freezer, 

its  thermal  expansion  coefficient  can  be  calculated  as: 

AL  nOxlO^m  , 

a  =  =  =  38xlO"6/°C 

LxAT    0.12604m  x  25°  C 


Table  6-1  Physical  properities  of  the  HCMM  components 


base  sphere 
ballscrew 
pin  follower 

triangle  frame 
cylinder  column 

strut  shoulder 

center  rod 

material 

stainless  steel 

cast  steel 

aluminum 

fiberglass 

density 
P 

(kg/m3) 

7750 

7700 

2750 

Not  Available 

thermal  conductivity 
k  @  300°K 
(W/m/  °C) 

17.5 

43 

190 

Not  Available 

specific  heat 
c 

(J/kg/°C) 

500 

440 

900 

Not  Available 

thermal  expansion 
coefficient 
a 

16xl0"6/°C 
(or8.9xl0"6/°F) 

14xlO"6/°C 
(or7.8xl0"6/°F) 

23.5xlO'6/°C 
(or  13xl06/°F) 

38xl0"6/°C 
(or  22.2xlO"6/°F) 

The  Biot  number  was  first  checked  for  verifying  the  lumped  system  model,  and  its 
values  for  each  thermally  sensitive  component  is  listed  in  Table  6-2.  It  has  to  be  noted 
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that  the  value  of  convection  coefficient  of  hc  is  dependent  on  the  type  of  fluid  and  the 
environmental  conditions,  such  as  the  flow  speed  of  the  fluid  over  the  surface  of  the  solid. 
In  general,  for  a  given  shop  floor  heating  or  cooling  system,  the  natural  convection 
coefficient  (hc)  is  between  10  and  25  W/m2/°C.  Here  we  assume  hc=15  W/m2/°C  for  the 
case  study  in  this  chapter.  Since  the  Biot  number  for  all  the  components  are  far  less  0.1, 
the  lumped  system  model  will  be  applied. 


Table  6-2  The  Biot  number  of  each  thermally  sensitive  component 


hc  (W/m2/°C) 

Ls  (m) 

k  (W/m/°C) 

Bi  =  hcxLs/k 

sphere 

15 

r/3=0.00847 

17.5 

0.0073 

pin  follower 

15 

r/2=0.00159 

17.5 

0.0016 

triangle  frame 

15 

w  /4=  0.0238 

43 

0.0083 

cylinder  column 

15 

r  12=  0.04286 

43 

0.0149 

For  a  lumped  system,  the  solution  for  transient  thermal  convection  can  be 
obtained  from  Equation  6-6.  Equation  6-6  tells  us  that  the  temperature  at  any  moment  in 
a  solid  component  is  a  function  of  its  initial  temperature,  its  thermal  time  constant,  and 
the  ambient  temperature.  The  thermal  time  constant  of  a  solid  component  will  determine 
the  quickness  of  the  system  in  reaching  equilibrium.  From  Figure  6-5(a)  we  may  observe 
that  a  component  with  small  time  constant  can  quickly  reach  equilibrium  with  ambient 
temperature.  Therefore,  a  small  thermal  time  constant  component  can  follow  the  ambient 
temperature  more  closely  than  a  component  with  large  time  constant  if  the  ambient  is 
under  a  dynamic  temperature  variation.  This  phenomena  will  be  demonstrated  later 
during  the  case  study.  Since  the  thermal  time  constant  is  such  an  important  parameter  for 
determining  the  thermal  response  of  a  solid  in  a  dynamic  environment,  followed  by  the 


90 


Biot  number,  the  time  constant  of  each  component  is  then  estimated  by  using  Equation  6- 


•  The  time  constant  for  reference  spheres: 

pVc     7750x%7rr3  x500 
T  .  „  =  =  :  =  2187.2  seconds  =  36.45  minutes 

sphere        ^  15x4^" 2 

•  The  time  constant  for  pin  followers: 

pVc     7750 xnr2Lx  500 
W  =  JJs  =  =  410.1  seconds  =  6.84  minutes 

•  The  time  constant  for  triangle  frame: 

pVc     7700  xw2Lx  440 
*triang,e  =        =  15x4wL  =  5378.5  seconds  =  89.6  minutes 

•  The  time  constant  for  column  cylinder: 

pVc     7700  x  (7Tr02  -  Ttr2  )L  x  440 
^column  =  T—r  =  i<x,o^.  t  =  20329  seconds  =  33.88  minutes 

Because  most  of  the  physical  properties  for  fiberglass  tube  are  not  available,  its 
time  constant  has  to  be  obtained  in  another  manner.  In  addition,  because  of  the 
uncertainty  existing  in  the  listed  physical  properties,  it  will  be  better  if  we  can  verify  the 
theoretic  thermal  time  constant  values  for  all  components.  In  order  to  obtain  the  time 
constant  for  fiberglass  tube  and  to  verify  the  theoretical  time  constant  value  for  other 
components,  several  thermal  experiments  have  been  performed. 
Time  Constant  Verification 

If  the  ambient  temperature  and  the  temperature  on  the  solid  component  can  be 
measured,  the  thermal  time  constant  of  the  solid  can  be  evaluated  by  using  Equation  6-7. 
Thus  we  performed  this  experiment  by  measuring  the  ambient  temperature  in  the 
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metrology  laboratory  where  the  HCMM  was  located.  At  the  same  time,  the  temperatures 
on  the  surfaces  of  the  solids  also  were  recorded.  The  thermistors  used  as  the  temperature 
sensors  have  a  measurement  resolution  of  0.0 1°C.  One  thermal  sensor  was  used  to 
measure  the  ambient  temperature.  The  surface  temperature  of  solid  components  are  taken 
by  averaging  the  three  thermal  sensors  which  were  placed  at  different  surface  locations. 
The  recorded  temperatures  for  different  components  are  plotted  in  Figure  6-7  to  Figure  6- 
9.  The  thermal  time  constant  for  each  solid  component  was  then  evaluated  by  using 
Equation  6-7.  From  Equation  6-7,  the  temperature  of  the  solid  was  estimated  by  plugging 
in  the  measured  ambient  temperature  (Tair)  and  guessing  a  time  constant  (x).  The  time 
step  (At)  in  Equation  6-7  is  the  time  interval  used  to  record  the  ambient  temperature 
which  is  0.5  minute  in  the  experiments.  If  the  assumed  time  constant  is  close  to  the  true 
value,  the  calculated  temperature  will  be  close  to  the  measured  temperature  in  the  solid. 
By  using  visual  comparison,  the  approximate  time  constants  are  38,  85,  and  35  minutes 
for  reference  sphere,  triangle  frame,  and  supporting  column  respectively.  Comparing  the 
experimental  time  constant  values  with  the  theoretical  values,  the  difference  are  within 
5%  .  We  then  have  the  confidence  to  use  the  values  of  the  theoretical  time  constant  for 
the  later  analysis. 

When  we  estimated  the  time  constant  for  fiberglass  tube,  it  was  assumed  that  the 
thermal  system  is  a  lumped  system.  Therefore,  the  method  mentioned  above  can  be 
applied.  The  time  constant  for  the  fiberglass  tube  was  estimated  to  be  6.8  minutes  as 
shown  in  Figure  6-10. 

Because  the  pin  follower  was  not  available  when  the  experiment  was  taken,  its 
time  constant  could  not  verified. 
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Estimated  Time  constant  for  Reference  Sphere  =  38minutes 

Temp  (°F) 
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Figure  6-7  Time  constant  estimation  for  Reference  Sphere 


Figure  6-8  Time  constant  estimation  for  triangle  frame 


93 


Estimated  Time  Constant  for  Column  :35  minutes 
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Figure  6-9  Time  constant  estimation  for  column 


Figure  6-10  Time  constant  estimation  for  fiberglass  tube 


94 

Thermal  Response  Under  Different  Thermal  Cycle 

Several  simulations  have  been  performed  to  help  us  understand  the  temperature 
relationship  between  the  ambient  and  the  thermally  sensitive  components  of  the  HCMM. 
In  our  simulations,  it  was  assumed  that  the  initial  temperature  in  each  component  is  the 
same  as  the  ambient.  The  ambient  temperature  then  varied  in  a  sinusoidal  manner  which 
is  written  as: 

D  2m 

Tair(t)  =  T0±-sm(—)  (6-9) 

where  T0  is  the  initial  ambient  temperature,  P  is  the  temperature  cycling  period, 
and  D  is  the  amplitude  of  temperature  variation.  If  the  ambient  temperature  variation 
amplitude  (D)  and  cycling  period  (P)  are  given,  the  temperature  in  each  component  at  any 
moment  can  be  calculated  by  using  Equation  6-7.  Figure  6-11  to  Figure  6-14  show  some 
of  the  simulation  results  in  which  the  ambient  is  under  different  cycling  frequency.  It  can 
be  observed  from  these  four  figures  that  the  faster  the  thermal  cycle  in  ambient 
temperature  the  smaller  the  temperature  change  in  solid  components.  This  phenomena 
gives  us  the  reason  why  in  a  good  metrology  laboratory  the  air  temperature  cycle 
frequency  is  often  very  fast  (for  instance  one  cycle  per  minute).  Another  conclusion 
which  can  be  made  for  the  thermally  sensitive  components  is  that  the  larger  the  time 
constant  the  smaller  the  temperature  variation  in  their  body. 
Thermal  Growth 

The  maximum  temperature  variation  in  a  solid  component  can  be  obtained  if  the  ambient 
temperature  frequency  and  its  variation  amplitude  (D)  are  given.  These  results  are 
tabulated  and  plotted  in  Table  6-3  and  Figure  6-15.  In  Table  6-3  the  values  given  are 
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Figure  6-1 1  Temperature  variation  for  5  cycle/hour 
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Figure  6-12  Temperature  variation  for  2  cycle/hour 


96 


Figure  6-13  Temperature  variation  for  1  cycle/hour 


Figure  6-14  Temperature  variation  for  0.5  cycle/hour 
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Table  6-3  The  maximum  temperature  change  in  components  for 
different  air  temperature  frequency 

Unit:  °F 


Temp.  freq. 
(cycle/hr) 

max(AT) 
(pin  follower, 
fiberglass  tube) 

max(AT) 
(column) 

max(AT) 
(sphere) 

max(AT) 
(triangle  frame) 

5 

0.3067 

0.0595 

0.0551 

0.0216 

2 

0.6115 

0.1561 

0.1447 

0.0567 

1.5 

0.7466 

0.2339 

0.2178 

0.0867 

1 

0.8256 

0.3072 

0.2867 

0.1171 

0.5 

0.9429 

0.5331 

0.5070 

0.2359 

0.25 

0.9846 

0.7668 

0.7451 

0.4352 

The  values  are  based  on  the  maximum  ambient  temperature  variation,  D=1°F 


Temp. 
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Figure  6-15  The  maximum  temperature  change  in  components  for 
different  air  temperature  frequency 
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normalized  against  a  value  of  D=1°F  .  The  number  found  in  the  table  then,  multiplied  by 
the  actual  value  of  maximum  ambient  temperature  variation,  D(°F ),  will  give  the 
maximum  temperature  for  the  respective  components.  If  we  assume  that  the  temperature 
change  in  components  exhibits  a  linear  relationship  with  time,  though  this  is  not  quite 
correct,  the  temperature  rate  of  change  in  each  component  can  be  obtained  by  taking  the 
maximum  temperature  change  in  the  component  divided  by  the  time  period  of  a  half 
temperature  cycle.  These  temperature  rates  of  change  are  listed  and  plotted  in  Table  6-4 
and  Figure  6-16. 

The  maximum  length  change  in  components  can  be  obtained  by  converting  the 
maximum  temperature  change  values  in  Table  6-3  by  using  the  thermal  expansion 
equation,  Equation  6-8  .  The  lengths  or  diameter  used  in  thermal  expansion  calculation 
for  fiberglass  tube,  pin  follower,  reference  sphere,  triangle  frame,  and  supporting  column 
are  33",  1",  2",  68.33"  and  33.46"  respectively.  The  thermal  expansion  coefficient  for 
each  component  can  be  obtained  from  Table  6-1.  By  solving  the  thermal  expansion 
equation,  we  have  the  maximum  length  changes  in  components,  and  their  values  are 
listed  and  plotted  in  Table  6-5  and  Figure  6-17.  Similarly  the  length  changing  rate  in 
each  component  can  be  obtained  by  converting  the  values  of  temperature  rate  of  change 
listed  in  Table  6-4.  After  performing  the  thermal  expansion  conversion,  the  length  rate  of 
change  is  tabulated  in  Table  6-6  and  plotted  in  Figure  6-18. 
Time  Period  Determination  for  Re-Calibration  or  Finishing  a  Part  Measurement 
The  values  of  the  maximum  length  change  and  the  length  growth  rate  provides  us  with 
valuable  information  for  determining  how  fast  we  have  to  finish  a  part  measurement  or 
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how  often  the  HCMM  has  to  be  re-calibrated  before  the  thermal  effect  exceeds  the 
allowable  value. 

For  instance,  the  HCMM  is  placed  in  the  metrology  laboratory  has  an  ambient 
temperature  variation  of  1 .4°F  and  the  frequency  of  2  cycles  per  hour.  If  the  maximum 
allowable  thermal  growth  in  length  for  the  triangle  frame  and  the  supporting  column  is 
43|X  inch  for  each,  we  want  to  understand  how  long  it  may  take  for  the  triangle  frame  or 
the  supporting  column  to  growth  more  than  43jll  inch,  such  that  the  re-calibration  process 
has  to  be  performed  or  a  part  measurement  has  to  be  finished.  In  this  case,  from  Table  6- 
5  we  found  the  maximum  length  change  in  the  triangle  frame  is  42.14(1  inch.  This  tells  us 
that  the  length  change  in  the  triangle  frame  during  the  whole  operation  period  will  remain 
in  the  allowable  range.  Therefore,  the  re-calibration  process  is  not  necessary  to  for  this 
component.  For  the  supporting  column,  its  maximum  length  change  will  be  56.84(1  inch 
which  is  more  than  the  allowable  thermal  growth  of  43\i  inch.  To  find  the  time  period  for 
the  supporting  column  to  growth  to  43jli  inch,  one  can  use  the  length  growth  rate  listed  in 
Table  6-6.  It  will  be  found  that  the  thermal  growth  rate  in  the  supporting  column  is  3.8 1  jLL 
inch/min  (i.e.  1.4°F  x2.72fx  inch/min/°F  ).  Therefore,  the  time  period  for  the  supporting 
column  to  grow  to  43(4,  inch  will  be  1 1.3  minutes.  A  part  measurement  operation 
therefore  has  to  be  finished  in  11.3  minute.  This  limited  operational  time  may  be  a 
problem  for  users  of  the  HCMM.  This  problem  can  be  improved  by  putting  a  thermal 
isolator  around  the  surface  of  the  supporting  column  which  will  reduce  the  temperature 
variation  range  inside  the  column,  or  by  filling  sand  into  the  hollowed  column  cylinder 
which  may  increase  the  thermal  time  constant. 
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Table  6-4  The  temperature  rate  of  change  in  each  component  for 
different  air  temperature  frequency 

Unit:  °F/min 


thermal  cycle 
(cycle/hr) 

AT 

(pin  follower 
fiberglass  tube) 

AT 
(column) 

AT 
(sphere) 

AT 

(triangle  frame) 

5 

0.05111 

0.00992 

0.00918 

0.00359 

2 

0.04077 

0.01041 

0.00965 

0.00378 

1.5 

0.03318 

0.01040 

0.00968 

0.00385 

1 

0.02752 

0.01024 

0.00956 

0.00390 

0.5 

0.01572 

0.00888 

0.00845 

0.00393 

0.25 

0.00821 

0.00639 

0.00621 

0.00363 

The  values  are  based  on  the  amplitude  of  air  temperature  variation,  D=1°F 


Figure  6-16  The  temperature  rate  of  change  in  each  component  for  different 
air  temperature  frequency 
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Table  6-5  The  maximum  length  change  in  components  for  different  air  temperature 
frequency 

 unit:u.  inch 


thermal  cycle 
(cycle/hr) 

max(AL) 
pin  follower 

max(AL) 
fiberglass  tube 

max(AL) 
column 

max(AL) 
sphere 

max(AL) 
triangle  frame 

5 

2.7 

231.64 

15.5 

1.0 

11.5 

2 

5.4 

425.98 

40.6 

2.6 

30.1 

1.5 

6.6 

520.09 

60.9 

3.8 

46.1 

1 

7.3 

575.13 

80.0 

5.0 

62.2 

0.5 

8.4 

656.86 

138.7 

9.0 

125.4 

0.25 

8.8 

685.92 

199.6 

13.2 

231.3 

The  values  are  based  on  the  amplitude  of  air  temperature  variation,  D=1°F 


cycle/hr 


Figure  6-17  The  maximum  thermal  growth  in  components  for 
different  thermal  cycle  with  air  temperature  variation  of  1°F 
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Table  6-6  The  length  growth  rate  in  components  for  different  thermal  cycle  rates 


Unit:  \i  inch/min 


thermal  cycle 
(cycle/hr) 

AL 

(pin  follower) 

AL 

(fiberglass  tube) 

AL 
(column) 

AL 
(sphere) 

AL 

(triangle  frame) 

5 

0.45 

35.61 

2.59 

0.16 

1.92 

2 

0.36 

28.40 

2.72 

0.18 

2.01 

1.5 

0.30 

23.12 

2.71 

0.18 

2.05 

1 

0.24 

19.17 

2.67 

0.18 

2.08 

0.5 

0.14 

10.95 

2.32 

0.16 

2.10 

0.25 

0.07 

5.72 

1.67 

0.12 

1.93 

The  values  are  based  on  the  amplitude  of  air  temperature  variation,  D=1°F 
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Figure  6-18  The  length  rate  of  change  in  components  for  different  thermal 
cycle  rates 
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The  maximum  Allowable  Calibration  Period 

In  the  opposite  way,  if  the  allowable  thermal  errors  in  the  strut  length  and  in 
center  rod  length  have  been  specified,  one  can  determine  the  maximum  allowable  time 
period  during  a  calibration  process  must  be  completed.  It  can  be  determined  by  choosing 
the  smaller  one  from  the  following  two  time  periods: 

allowable  AL„ 


Times,ru, 


strut 


Strut  length  changing  rate 


allowable  ALcenlerroJ 

Tlme  center  = 


Center  rod  length  changing  rate 
For  the  thermal  environment  mentioned  above,  if  the  maximum  allowable  thermal 
error  for  the  unsensed  region  in  strut  is  4fx  inch  and  for  the  center  rod  length  error  is  20|i 
inch.  Then  the  maximum  allowable  calibration  time  will  be  the  shorter  one  from  the 
following  two  values 

Timestrut=  4n'7[1.4°Fx(0.36x2+0.18)  u."/°F.min  ]  =  3.17  minutes 
Time center  Wd  =  20|jy[  1 .4°Fx(28.40+0. 1 8)  n'7°F.min  ]  =  30  seconds 
In  this  case  the  maximum  allowable  time  period  for  completing  a  calibration  procedure  is 
30  seconds.  This  time  period  could  be  too  short  for  us  to  finish  a  calibration  procedure. 
To  prolong  the  time  period,  we  suggest  changing  the  fiberglass  tube  into  a  material  that 
has  the  smallest  possible  thermal  expansion  coefficient  such  as  Zerodur. 

6.3  The  Effect  of  the  Heat  Generated  in  Actuators 

The  actuators  are  expected  to  have  a  performance  of  the  maximum  acceleration  of 
64.4ft/sec2  (or  2G's),  a  maximum  linear  speed  of  20  in/sec  (or  0.5m/sec),  and  a  maximum 
expected  load,  according  Mevoli  [4],  of  1501b  (or  667N).  It  is  known  that  a  motor  will 
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generate  heat  during  operation  due  to  the  friction  and  electronic  resistance.  Thus,  through 
heat  conduction,  the  heat  in  the  actuator  will  be  transmitted  to  the  ballscrew  and  the 
shoulder.  Let  us  consider  an  extreme  HCMM  operation  case:  A  constant  strut  extension 
speed  of  20  in/sec  (0.5m/sec)  was  maintained  over  the  travel  range,  and  the  actuator  was 
subjected  to  an  axial  load  of  1501b  (667N)  during  the  entire  process.  Before  the  actuator 
started  to  function,  the  whole  strut  system  had  the  same  temperature  as  the  ambient.  Here 
the  ambient  temperature  is  assumed  to  be  constant  during  the  whole  operation  process. 
For  this  analysis,  we  assume  that  the  motor  loses  33%  of  its  power  into  heat  and  only  uses 
its  67%  of  the  power  to  generate  the  strut  motion.  In  this  section  we  will  investigate  that 
how  heat  inside  the  motor  will  be  transmitted  to  the  pin  followers,  and  result  in  their 
length  changing. 

A  power  of  333.5  Watts  (N.m/sec),  which  is  the  axial  load  multiplied  by  the  strut 
speed,  has  to  be  consumed  to  complete  the  required  motion.  From  our  assumption  a  heat 
rate  of  167.5  Watts  will  be  generated  in  the  motor.  To  evaluate  the  temperature  (or  the 
heat)  at  the  end  of  ballscrew  or  strut,  a  finite  element  method  was  used.  The  finite 
element  (FE)  model  is  made  as  shown  in  Figure  6-19.  Because  of  the  symmetric 
geometry  in  the  strut,  we  only  consider  half  of  our  investigating  region.  Since  there  is  no 
temperature  gradient  along  the  symmetry  line,  the  boundary  condition  along  this 
symmetry  line  is  treated  as  if  it  is  insulated.  Thermal  convection  is  considered  between 
the  strut  surface  and  the  ambient,  and  its  convection  coefficient  (hc)  is  assumed  to  be  15 
W/m  /°C.  The  transient  analysis,  then,  was  performed  from  time  t=0  to  t=20  minutes  by 
using  a  commercially  available  FE  software  COSMOS/M. 
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The  resulting  temperature  distributions  for  time=l  minute  and  time=20  minutes 
are  plotted  in  Figure  6-20  and  Figure  6-21.  It  can  be  observed  that  the  temperature  inside 
the  motor  is  raised  nearly  by  12  °C  after  20  minute's  operation.  The  temperature  history 
at  each  end  of  the  ballscrew  and  the  shoulder  are  shown  in  Figure  6-22  and  Figure  6-23 
respectively.  Because  the  thermal  convection  occurred  between  the  air  and  the  surface  of 
the  strut,  it  can  be  expected  that  most  of  the  transmitted  heat  from  the  motor  will  be 
dissipated  into  air  before  it  reaches  the  end  of  the  ballscrew  or  shoulder. 

It  can  be  observed  that  there  is  almost  no  temperature  raise  within  the  first  four 
minutes.  This  implies  that  the  motor's  heat  will  not  affect  the  pin  followers  length 
change  during  the  calibration  process  if  this  process  can  be  completed  in  4  minutes. 

6.4  Conclusion 

The  thermal  environment  is  such  a  complicated  issue  that  it  is  difficult  to  predict 
what  the  ambient  temperature  would  be.  In  this  chapter  we  provide  a  methodology  for 
evaluating  the  allowable  time  period  before  the  HCMM  needs  to  be  re-calibrated  or  a  part 
measurement  has  to  be  finished,  and  how  long  a  calibration  process  can  last.  This 
method  is  based  on  knowledge  of  the  environmental  conditions.  In  order  to  have  this 
information,  it  is  required  to  monitor  the  ambient  temperature  in  the  working  area. 

We  also  suggest  to  change  some  design  parameters  in  the  supporting  column  and 
the  center  rod  so  that  the  machine  operation  time  and  the  calibration  processing  time  can 
be  prolonged. 
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Figure  6-19  The  finite  element  model  for  a  strut 
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Figure  6-20  The  temperature  distribution  in  strut  after  the 
motor  actuates  for  1  minute 
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Figure  6-21  The  temperature  distribution  in  the  strut  when 
motor  actuates  for  20  minutes. 
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Figure  6-22  The  temperature  history  at  the  end  of  ballscrew 


Figure  6-23  The  temperature  history  at  the  end  of  the  shoulder 
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Figure  6-24  The  temperature  distribution  along  the  ballscrew  at  t=20 
minutes 


CHAPTER  7 

HCMM  MEASUREMENT  UNCERTAINTY  ANALYSIS 

In  this  chapter  we  will  identify  the  error  sources  that  will  affect  the  HCMM 
measurement  accuracy,  and  investigate  how  the  error  sources  propagate  through  the 
mathematical  model  of  the  HCMM  system  to  assess  how  they  impact  the  probe 
coordinate  calculations  and  cause  errors  in  the  measured  position  of  the  probe. 

Figure  7-1  sketches  the  major  sources  that  affect  the  determination  of  the  probe 
coordinates.  Three  major  error  sources  will  be  discussed  in  this  chapter:  linear 
displacement  measurement  error,  calibration  error,  and  thermal  error.  Other  error  sources 
such  as  the  probe  alignment  error  and  the  errors  due  to  structural  vibration  will  not  be 
included  in  our  investigation,  since  they  are  beyond  the  scope  of  this  research. 

7.1  Error  Sources  In  Strut  Displacement  Measurement 

Four  major  error  sources  are  considered  to  be  the  contributors  for  the  strut  linear 
displacement  measurement  error: 

1 .  Laser  interferometer  error 

2.  Spherical  joint  imperfection 

3.  Reference  sphere  imperfection 

4.  Thermal  growth  in  the  pin  followers  and  spheres 

7.1.1  Laser  Interferometer  Error 

Three  major  error  contributors  for  the  laser  interferometer  measurement 
uncertainty  are  identified  as  follows, 
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imperfection 


reference 
sphere 
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Figure  7-1  The  error  sources  that  affect  the  probe  position  determination 
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1 .  Optical  component  rotation 

2.  Refractive  index  error  and  environmental  effects 

3.  Laser  stability  and  intrinsic  factor 

Optical  Component  Rotation 

This  error  comes  from  misalignment  of  the  laser  beam  with  the  extension 
direction  as  a  result  of  the  rotation  of  the  emitter  and  retroreflector  at  either  end  of  the 
strut.  The  rotation  of  the  optics  is  caused  by  compliance  in  the  bearings  and  flexure  (sag) 
of  the  strut  tubes  due  to  gravity.  Figure  7-1  illustrates  the  optical  component  rotation 
error.  According  the  analysis  of  Mevoli  [4],  the  strut  designer,  this  error  is  between 
0-2.4H  inch  (0.06245^im). 


Non-Deflected  Strut 


23 


Actual   Intended  target 

Target  " 


Figure  7-2  Optical  component  rotation  error 


Refractive  Index  Error  and  Environmental  Effects 

Temperature,  pressure,  humidity,  and  gas  composition  affect  the  refractive  index 
of  light  in  the  measurement  area.  Changes  in  the  refractive  index  causes  a  cumulative 
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error  over  the  path  of  measurement  beam.  Because  the  vacuum  system  will  be  installed 
on  the  HCMM  the  strut  will  be  evacuated.  The  environmental  effect  will  be  dramatically 
reduced.  We  estimate  that  the  uncertainty  caused  by  the  refractive  index  error  and 
environmental  will  be  below  the  range  of  ±luin. 
Laser  Stability  and  Intrinsic  Factors 

The  stability  of  the  wave  length,  the  electrical  noise  error  in  the  phase  detector, 
and  the  optical  components  (the  surface  finish  effects  and  coating  property  variations  on 
the  polarization  beamsplitter)  are  inherent  in  all  heterodyne  interferometry  systems.  The 
magnitudes  however  are  unique  to  particular  laser  systems.  For  the  Zygo  Axiom  2/20 
laser  system,  these  errors  are  in  the  range  of  ±0.63|i  inch. 
7.2.2  Spherical  Joint  Imperfection 

The  strut  length  is  defined  as  the  length  between  the  base  sphere  center  to  the 
moving  sphere  center.  If  the  slide  motion  is  not  parallel  to  the  line  between  the  two 
spheres,  the  measured  strut  displacement  will  not  be  correct.  This  error  is  caused  by  the 
imperfection  of  the  spherical  joint.  Because  of  the  manufacturing  tolerances  on  joint 
components,  installation  error,  and  the  joint  deflection  caused  by  the  static  or  dynamic 
load,  the  spherical  joint  inevitably  will  be  off  center  during  the  spherical  motion.  Figure 
7-3  shows  the  linear  displacement  measurement  error,  R-R',  caused  by  the  joint  offset  of 
(b,h)  in  the  X  and  Y  directions  from  the  perfect  joint  point.  To  calculate  the  measurement 
error,  four  assumptions  are  made:  (1)  the  pin  follower  is  aligned  with  the  laser  beam,  (2) 
the  end  face  of  the  pin  follower  is  always  tangent  to  the  reference  sphere,  (3)  reference 
sphere  is  perfect  sphere,  (4)  the  end  face  of  pin  follower  is  perpendicular  to  its  axial 
direction. 
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Joint  offset  in  X  direction  (b) 


R   -v-  rj 


Figure  7-3  Spherical  joint  imperfection 


The  notations  show  in  Figure  7-3  are  defined  as: 
O:  The  center  of  base  sphere. 

R:  The  distance  between  base  sphere  center  to  retroreflector  (without  joint  offset). 
R':  The  distance  between  base  sphere  center  to  retroreflector  (with  joint  offset), 
rj:  The  spherical  joint  radius 

v:  The  distance  from  spherical  joint  to  retroreflector. 
Pc:  The  contact  point  of  the  pin  follower  on  the  sphere  (without  joint  offset). 
Pc':  The  contact  point  of  the  pin  follower  on  the  sphere  (with  joint  offset).  Its 
coordinate  is  (xc,yc). 

Pe:  the  point  that  the  extended  laser  beam  hits  on  the  end  face  of  the  pin  follower 

(without  joint  offset). 
Pe' :  the  point  that  the  extended  laser  beam  hits  on  the  end  face  of  the  pin  follower 

(with  joint  offset).  Its  coordinate  is  (xuyi) 
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b:  The  joint  offset  in  the  x  direction 
h:  The  joint  offset  in  the  y  direction, 
r:   The  radius  of  the  sphere. 

To  find  the  locations  of  Pe'  and  Pc',  the  following  four  geometrical  constraints  are  used 

1.   Line  equation  for  OPe'  : 

„         h  y, 

R  +  v-b  jc, 


2.  Line  equation  for  Pe'Pc' 


tan(0  +  9O°)  =  - — —  (7-2) 


*,  -  xc 


3.   Radius  of  sphere: 


r  =  ^(R  +  v  +  rj-xc)2  +  y2  (7-3) 

4.   Line  Pe'Pc'  tangents  to  sphere: 

( R  +  v  +  rj  -  x}  f  +  y2  =  {(*,  -  xc  f  +  ( y,  -  yc )2}  +  r2  (7-4) 

For  the  case  R=41",  rj=1.5",  v=0.5",  r=l",  b=0.03",  h=0.03",  the  coordinates  of 
Pc'  and  Pe'  can  be  found  by  using  Equations  7-1  to  7-4,  and  their  locations  are 
(42.00000026,  -0.0007234147)  and  (41.99997779,0.0303394976)  respectively.  Thus  the 
linear  displacement  measurement  error  (R-R')  caused  by  the  joint  offset  of  (0.03",  0.03") 
will  be  1 1.26|i  inch.  At  different  joint  offset  magnitudes  the  linear  displacement 
measurement  errors  are  calculated  and  listed  in  Table  7-1.  From  Table  7-1  it  can  be 
observed  that  the  major  measurement  error  is  caused  by  the  joint  offset  in  the  Y-direction 
(the  direction  normal  to  the  line  which  passes  through  the  two  sphere  centers) .  The 
relationship  between  the  joint  offset  in  Y-direction  and  the  measurement  error  is  plotted 
in  Figure  7-4. 

For  the  spherical  joints  we  built,  the  measured  maximum  joint  offset  in  the  radial 
direction  of  the  strut  is  0.002  inch  according  to  the  experiments  of  Yeager  [47].  Thus  we 
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can  very  safely  assume  that  the  measurement  uncertainty  caused  by  the  joint  offset  will  be 
less  than  0.62|iin  (0.31|iin  for  each  joint). 


Table  7-1  The  linear  displacement  error  caused  by  joint  offset 


joint  offset  (b,h) 

linear  displacement 

(inch) 

measure  error  (inch) 

(0.005,  0.005) 

3.10E-07 

(    0  ,0.005) 

3.10E-07 

(0.01,0.01) 

1.25E-06 

(  0  ,0.01) 

1.25E-06 

(0.03,  0.03) 

1.13E-05 

(  0  ,0.03) 

1.124E-05 

(0.05,  0.05) 

3.13E-05 

(  0,0.05) 

3.12E-05 

(0.1,0.1) 

1.25E-04 

(  0,0.1) 

1.2485E-04 

140 


0        0.02      0.04      0.06      0.08  0.1 
joint  offset  in  the  normal  direction  (inch) 


Figure  7-4  The  measurement  error  caused  by  joint  offset 


117 

7.1.3  Reference  Sphere  Imperfection 

The  reference  sphere  is  similar  to  a  master  gage  in  that  it  provides  a  constant 
reference  length  (i.e.  radius)  from  its  surface  to  center.  The  imperfections  of  form 
(sphericity)  of  this  sphere  will  cause  uncertainty  on  the  reference  length.  For  the  HCMM 
the  Grade  8  spheres  with  sphericity  of  8  (xin  were  used.  The  sphericity  is  defined  as  the 
radius  difference  between  the  maximum  inscribed  and  the  minimum  circumscribed 
spheres  (as  shown  in  Figure  7-5).  The  maximum  reference  length  deviation  could  be  8 
\iin.  However,  the  region  used  as  a  reference  is  in  a  area  of  60°  angle,  we  assume  the 
reference  length  deviation  in  this  contact  region  will  be  within  ±4|xin. 


Figure  7-5  The  sphericity  of  a  sphere 

7. 1 .4  Thermal  Error 

In  a  strut,  the  laser  interferometer  system  will  detect  a  length  change  from  its 
emitter  to  the  retroreflector,  as  shown  in  Figure  6-1.  However,  the  thermal  growth  in  the 
components  (pin  followers  and  spheres)  which  is  not  covered  by  the  laser  path  will  not  be 
detected.  The  unmeasured  thermal  growth,  then,  will  be  treated  as  measurement  error  in 
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the  strut  length.  The  environment  in  which  the  HCMM  is  located  has  an  ambient 
temperature  frequency  of  2  cycles  per  hour  and  a  maximum  variation  amplitude  of  1 .4°F. 
From  Table  6-6,  the  thermal  growth  for  the  unsensed  region  will  be  within  ±1.26  uJn/min 
(  1.4x(2x0.36+0.18). 


Table  7-2  error  sources  encountered  for  strut  linear  displacement  measurement 


error  sources 

range  of  error 
(uin) 

1 

Laser  interferometer 

Optical  component  rotation 

0-2.4 

2 

Deadpath  and  environmental  effects 

-1-+1 

3 

Laser  stability  and  intrinsic  factor 

-0.63-+0.63 

4 

Spherical  joint  imperfection 

-0.62-0.62 

5 

Reference  sphere  imperfection 

-8~+8 

6 

Thermal  growth  in  pin  followers  and  spheres  in  1  minute 

-1. 26-+ 1.26 

Error  summary 

-11.51-13.91 

1.2  The  Combined  Standard  Uncertainty 

Suppose  a  measurand  of  interest  F  is  related  to  n  independent  quantities  xi ,  i 

=  1 ,. . . ,n.  The  measurand  of  interest  F  can  be  written  as 

F  =  /(*,,...*„)  (7-5) 
To  determine  F,  the  quantities  x,  may  be  directly  measured  or  simply  estimated;  in 

either  case  our  state  of  knowledge  of  these  quantities  is  described  by  a  set  of  probability 

distributions  with  associated  variances  (J x  .  The  positive  square  root  of  these  variances 

are  called  "standard  uncertainty"  of  the  quantities  x,.  Using  the  Taylor  expansion,  the 
measurand  of  interest  F  can  be  rewritten  as 
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v-  df 

F  =  F0  +2^  3~  Ac,.  +  higher  order 


terms 


(7-6) 


where  the  F0  is  the  nominal  value  of  F,  which  is  obtained  by  plugging  the  nominal  values 
of  Xi  into  Equation  7-5.  If  the  nominal  values  of  x,  are  close  enough  to  true  values,  the 
higher  order  terms  can  be  ignored.  Thus,  the  error  associated  with  the  F0  can  be  rewritten 


as 


df 
,=t  dx, 


(7-7) 


Since  the  error  of  each  quantity,  Ac, ,  is  a  probability  distribution,  the  measured 
error  of  interest  F  must  be  a  probability  distributed  quantity.  In  order  to  evaluate  the 
quantity  of  AF,  the  statistical  method  has  to  be  used.  The  expected  value  for  the  quantity 
of  (AF)  can  be  written  as 


E[(AF)2]  =  E 


(tt  df 
y  — —  Ax 


(7-8) 


or 


m 

E[(AF)2]  =  ^ 


KdxtJ 


■E 


m-l  m 
i=l  j=i+l 


df  df 

dxt  dxj  j 


■E 


[Ax^j] 


(7-9) 


where  the  expectations  of  £[(AF)2]  and  £[Ac,2]  are  the  variances  of  the  quantities.  The 
expectation  of  £[Ax,.Ax;]  is  the  covariance  of  the  quantities.  Thus  the  "combined 
standard  uncertainty"  of  F  is  given  by  taking  the  positive  square  root  of 


»c(^)  =  Jl 


df 


m-l  m 


rd£d£^ 

ydxi  dx^j 


Cov(Ax,Axy.  \ 


(7-10) 
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For  the  HCMM,  if  the  standard  uncertainty  for  each  of  the  kinematic  parameters  is 
available,  the  combined  standard  uncertainties  for  upper  and  lower  sphere  positions  can 
be  estimated.  For  instance,  the  upper  sphere  at  position  i  in  the  X-direction,  Xup(i),  can  be 
written  as 


/2(U)-/2(3,/)  +  r2 


(7-11) 


The  combined  standard  uncertainty  for  the  coordinate  of  Xup(i),  then  can  be 
obtained  by  using  Equation  7-10.  It  can  be  written  as 


"c(*up): 


f(dX  V 

up 


<J,  + 


dx 


up 


dx 


.  1/2 


up 


dr  J 


2^^cov(A/,  A/3 )  +  2^^cov(A/,  Alr ) 


dlx     d  l3 

dl3    dlr  3  ' 


dlx  dlr 


V 


(7-12) 


or 


u  J  a"  + 


fi  ^ 


o,  + 


2r2 


+ 


+  2 


t,t3 


cov(A/,A/3)  +  2 


^  J 


-(V-v)  ,  i 


2r' 


+  — 
2 


cov(A/,  Ar) 


+  2 


(A  -V)  ,  1 
2r2  2 


cov(A/3Ar) 


1/2 


(7-13) 


From  Equation  7-13,  we  can  observe  that  the  sphere  position  uncertainty  is 
dependent  on  the  sphere's  location.  At  different  location,  the  strut  lengths  will  vary,  thus 
the  combined  standard  uncertainty  of  the  spheres  will  be  different. 


121 

To  determine  the  probe  position,  first  the  probe  orientation  to  the  center  rod  and 
the  offset  length  between  the  probe  tip  (stylus)  and  the  lower  sphere  center  has  to  be 
obtained.  In  a  conventional  CMM,  a  probe  verification  procedure  has  to  be  performed 
before  the  probe  can  be  used  to  make  measurement.  For  the  HCMM  this  procedure  also 
is  necessary  prior  to  the  measurement  operation.  Since  the  probe  verification  is  beyond 
the  scope  of  our  research,  here  we  assume  that  the  probe  is  aligned  and  parallel  with  the 
center  rod,  and  the  distance  between  its  stylus  and  lower  sphere  center  is  known  as  Lp. 


Upper  sphere 


Lc 


Lower  sphere 


Figure  7-6  The  center  rod  and  probe 


From  Figure  7-6,  the  probe  position  can  be  expressed  by  the  positional  vector, 


P      =  P    +  T  ■ 

1  probe       'low    '  ^p 


Lc 


(7-14) 


where  Lp  is  the  length  from  the  lower  sphere  center  to  the  probe  tip  (stylus).  Lc  is  the 
center  rod  length.  Therefore,  the  X-coordinate  of  the  probe  can  be  written  as 
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X  probe  ~  Xlow  +  ^p 


Xlow  Xup 


Lc 


1+- 


Lc 


Uow 


Lc 


up 


(7-15) 


Refer  to  Equation  7-7,  the  Axprobe  can  be  expressed  as 


Ax 


probe 


1  + 

Lc 


A/; 


A  L'fdx 


-— I 


up 


v<9/,. 


Lc 


2(xtow-x„p)ALc  (7-16) 


where  AI7,  Als,  and  A/9  equal  to  Ar,  Ab,  and  A/i  respectively.  Considering  the  last  term  in 
Equation  7-1,  if  the  center  rod  remains  parallel  to  the  Z  axis  (or  normal  to  the  granite 
table),  the  coordinates  of  Xiow  and  Xup  will  be  the  same.  Therefore,  the  last  term  can  be 
eliminated.  Equation  7-12  can  be  rewritten  as 


Ax 


probe 


(     1  \  » 

J  w 


1  +  -^- 
Lc 


'■low 


A/, 


L„  * 


 p_ 

Lc 


(d. 


up 


A/. 


(7-17) 


Thus,  the  variance  of  (xprobe  j  can  be  expressed  as, 


^7  A2 


Uc(X,ow)  + 


yLCj 


-2 


1  +  -^ 
Lc 


Lc 


,=1  .,=1 


Wow  «P 


(7-18) 


Therefore,  the  combined  standard  uncertainty  for  the  probe  position  in  the  X  direction  can 
be  obtained  by  taking  the  positive  square  root  of  Equation  7-18. 


UdX  probe) 


L  ^ 


Lc 


-2 


(1  \ 


UMlow)  + 


yLCj 


uSxuP) 


-*-YE 
Lc 


(dx.  d 


II 


XlowOXup  M 


v  dl,  dlj 


1/2 


(7-19) 


The  probe  uncertainties  in  the  Y  and  Z-directions  can  be  obtained  by  replacing  the 
character  "X"  in  Equation  7-19  by  "Y"  or  "Z".  For  our  current  HCMM  design,  the  probe 
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length,  Lp,  is  5  inches  and  the  center  rod  length  is  35  inches.  Therefore,  from  Equation  7- 
19,  we  can  observe  that  the  lower  sphere  position  uncertainty  will  be  the  major 
contributor  on  the  probe  position  uncertainty. 

7.3  Probe  Positional  Uncertainty  Due  to  Calibration  Error 

We  performed  two  thousand  calibration  simulations  by  using  the  optimal  pose  set 
(30  poses)  as  listed  in  Appendix  B.  In  each  simulation,  uniformly  distributed  noise  (- 
1 1.51 -+13.91  |iin)  was  superimposed  on  the  exact  strut  linear  displacement  to  simulate 
the  measurement  uncertainty.  From  the  2000  calibration  simulation  results,  we  plot  the 
histograms  of  the  identified  kinematic  parameter  errors  in  Figure  7-7.  It  can  be  observed 
that  the  identified  parameters  are  associated  with  an  unbiased  normal  distribution  error. 
To  represent  the  standard  uncertainty  on  each  of  the  kinematic  parameters,  and  to 
understand  the  covariance  between  the  parameters,  we  performed  the  descriptive 
statistical  data  analysis  by  using  the  calibration  simulation  results.  The  variance  and 
covariance  of  the  kinematic  parameters  are  listed  in  Table  7-3. 


Table  7-3  The  variance  of  identified  parameters  (m=30,  optimal  pose  set) 


al 

a2 

a3 

a4 

a5 

a6 

r 

b 

h 

al 

2.07E-08 

2.06E-08 

-8.10E-09 

-6.75E-09 

-1.23E-09 

-2.31E-09 

1.24E-08 

3.65E-08 

4.58E-10 

a2 

2.06E-08 

2.17E-08 

-7.11E-09 

-6.45E-09 

1.59E-09 

8.42E-11 

1.30E-08 

3.65E-08 

3.58E-09 

a3 

-8.10E-09 

-7.11E-09 

4.48E-08 

4.13E-08 

-1.20E-08 

-9.88E-09 

3.43E-08 

-3.74E-08 

-5.54E-09 

a4 

-6.75E-09 

-6.45E-09 

4.13E-08 

3.94E-08 

-9.30E-09 

-8.21E-09 

3.25E-08 

-3.43E-08 

-3.51E-09 

a5 

-1.23E-09 

1.59E-09 

-1.20E-08 

-9.30E-09 

5.90E-08 

5.37E-08 

-1.69E-08 

1.85E-09 

6.08E-08 

a6 

-2.31E-09 

8.42E-11 

-9.88E-09 

-8.21E-09 

5.37E-08 

4.99E-08 

-1.61E-08 

-6.66E-10 

5.58E-08 

r 

1.24E-08 

1.30E-08 

3.43E-08 

3.25E-08 

-1.69E-08 

-1.61E-08 

4.50E-08 

4.36E-10 

-9.87E-09 

b 

3.65E-08 

3.65E-08 

-3.74E-08 

-3.43E-08 

1.85E-09 

-6.66E-10 

4.36E-10 

7.86E-08 

7.11E-10 

h 

4.58E-10 

3.58E-09 

-5.54E-09 

-3.51E-09 

6.O8E-08 

5.58E-08 

-9.87E-09 

7.11E-10 

6.45E-08 
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Figure  7-7  The  histograms  of  identified  parameter  errors  (m=30) 
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The  variance  and  covariance  in  Table  7-3  provide  us  the  necessary  information  for 
calculating  the  combined  standard  uncertainty  for  the  probe  position.  As  we  know,  the 
combined  standard  uncertainty  will  vary  from  configuration  to  configuration.  When  the 
center  rod  remains  in  parallel  to  the  Z-axis,  at  each  configuration  the  combined  standard 
uncertainty  for  the  probe  can  be  obtained  by  using  Equation  7-19.  We  searched  the  entire 
HCMM's  work  volume  as  shown  in  Figure  7-8,  and  found  the  maximum,  minimum,  and 
the  mean  positional  uncertainties  for  the  probe,  as  listed  in  Table  7-4.  The  combined 
standard  uncertainties  are  plotted  in  Figure  7-9  where  the  probe  remained  in  the  plane  of 
Z=-25. 


Z  axis 


Y  Axis 


Figure  7-8  The  HCMM  work  volume 


Table  7-4  The  probe  positional  uncertainty  in  the  work  volume  due  to  the  calibration 
error 


m=30 

m=50 

uc(Probe_x) 

uc(Probe_y) 

uc(Probe_z) 

uc(Probe_x) 

uc(Probe_y) 

uc(Probe_z) 

mean 

166.4  uin 

59.2  uin 

33.9  uin 

118.1  uin 

46.8  uin 

25.1  ^in 

maximum 

190.3  uin 

69.5  |ain 

45.0  uin 

139.2  uin 

63.9  uin 

34.8  uin 

minimum 

147.4  uin 

54.8  uin 

27.4  uin 

101.7  uin 

39.2  uin 

21.3  uin 
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Probe  position  Uncertainty  ~Uc(Probe_x) 


Probe  position  Uncertainty  -Uc(Probe_y) 


Probe  position  Uncertainty  --Uc(Probe_z) 


Figure  7-9  The  combined  standard  uncertainty  for 
probe  at  the  level  of  Z=-25 
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Using  the  optimal  pose  set  with  50  configurations,  we  repeated  the  above 
mentioned  procedure,  the  probe's  positional  uncertainties  in  the  work  volume  were 
obtained  and  listed  in  Table  7-4. 

It  can  be  observed  from  Table  7-4  that  the  Z-direction  has  the  smallest  positional 
uncertainty  and  the  X-direction  has  the  largest  positional  uncertainty.  This  results  agree 
with  the  results  that  we  have  discussed  in  Section  5-4.  For  detailed  explanation  on  this 
phenomena  one  can  refer  back  to  Section  5-4. 

7.4  Probe  Position  Uncertainty  Due  to  Thermal  Error  and  Strut  Displacement 

Measurement  Error 

In  the  metrology  laboratory  where  the  HCMM  is  placed,  the  environment  has  an 
ambient  temperature  variation  of  1 .4°F  and  the  frequency  of  2  cycles  per  hour  (30 
minutes  a  cycle).  The  maximum  thermal  growth  of  the  triangle  frame  is  42.21  jxin.  Thus 
the  HCMM's  kinematic  parameter  "r"  will  be  subjected  to  a  range  of  error  within 
±42.2  luin.  For  the  default  values,  the  parameters  "ft"  and  "h"  will  be  1/2  and  V3 12 
times  of  the  value  of  "r."  For  a  uniformly  distributed  quantity  within  a  range  of  ±D,  the 
standard  uncertainty  will  be  D/V3  (see  Taylor  [42]).  Therefore,  the  variance  for  r,  ft,  and 
h  will  be  (42.21  x  10"6  )2  /  3  ,  [  42.21  x  10"6  x  - 1  /  3 ,  and 

422lxi0"6x—   / 3  respectively, 
k  2)  {  2) 

Since  the  variations  in  r,  b,  and  h  are  dependent  on  the  ambient  environment,  these  three 
parameters  must  be  completely  correlated  with  each  other.  The  covariance  between  the 


(r,  ft),  (r,  h),  and  (ft,  h)  then  can  be  written  as 


D  D 

x- 


S  2V3, 


D_  Sd) 
V3  2V3 


,  and 


(_D_xSp) 
2V3  X  2V3 


respectively,  where  D  is  the  value  of  42.21xlO"6  The  strut  linear  displacement 


128 

measurement  error  has  been  discussed  in  Section  7- 1 .  Here  we  will  use  the  results  in 
Table  7-2,  except  the  value  for  thermal  error  in  the  pin  followers  and  the  spheres  is  18.9 
pin  which  is  the  maximum  undetected  thermal  growth  in  the  strut.  The  error  range  in  the 
strut  displacement  measurement  then  will  be  within  ±31.55(J.in.  Thus  the  standard 
uncertainty  for  each  strut  will  be  31.55xlO~6/V3  .  Therefore,  the  variance  for  each  strut 
will  be  3.32xl0"10.  Since  there  are  many  error  sources  which  contributed  to  the  strut 
displacement  measurement,  it  will  be  very  difficult  to  determine  the  correlation  between 
struts.  Though  we  know  the  correlation  exists,  the  covariance  between  struts  ,for 
simplicity,  was  assumed  to  be  zero  which  means  the  error  in  each  strut  is  independent 
from  the  other.  The  zero  covariance  assumption  between  the  base  sphere  length  and  the 
strut  length  was  made  for  the  same  reason.  The  variance  and  covariance  between  the 
HCMM's  kinematic  parameters  are  then  listed  in  Table  7-5. 


Table  7-5  The  variance  and  covariance  between  kinematic  parameters 


li 

12 

1. 

u 

u 

16 

r 

b 

h 

1. 

3.32E-10 

0 

0 

0 

0 

0 

0 

0 

0 

h 

0 

3.32E-10 

0 

0 

0 

0 

0 

0 

0 

h 

0 

0 

3.32E-10 

0 

0 

0 

0 

0 

0 

u 

0 

0 

0 

3.32E-10 

0 

0 

0 

0 

0 

lj 

0 

0 

0 

0 

3.32E-10 

0 

0 

0 

0 

k 

0 

0 

0 

0 

0 

3.32E-10 

0 

0 

0 

r 

0 

0 

0 

0 

0 

0 

5.939E-10 

2.969E-10 

5.143E-10 

b 

0 

0 

0 

0 

0 

0 

2.969E-10 

1.485E-10 

2.572E-10 

h 

0 

0 

0 

0 

0 

0 

5.143E-10 

2.572E-10 

4.454E-10 

From  Table  7-5  we  have  the  necessary  information  for  determining  the  probe's 
positional  uncertainty.  By  using  Equation  7-19,  we  searched  the  HCMM's  work  volume, 
and  the  maximum,  minimum,  and  the  mean  positional  uncertainties  for  the  probe  were 
listed  in  Table  7-6.  We  observed  that  the  maximum  positional  uncertainty  is  in  the  Z- 
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direction.  Unlike  the  results  caused  by  the  calibration  error  as  shown  in  Table  7-4,  the 
kinematic  parameters  in  this  case  do  not  provide  a  special  relationship  to  minimize  the 
deviation  for  the  center  rod  length.  The  Z-coordinate  then  accumulates  the  uncertainties 
from  X-  and  Y-coordinate,  thus  in  the  Z-direction  will  have  the  largest  positional 
uncertainty. 


Table  7-6  The  combined  standard  uncertainties  for  the  probe  position  due  to  the  thermal 
error  and  strut  linear  displacement  measurement  error 


uc(Probe_x) 

uc(Probe_y) 

uc(Probe_z) 

mean 

23.3  uin 

20.3  uin 

62.9  uin 

maximum 

26.5  uin 

24.0  Uin 

99.6  uin 

minimum 

20.6  uin 

17.1  uin 

36.6  uin 

In  Figure  7-10  we  observed  that  the  probe's  positional  uncertainty  in  X- 
coordinate,  Uc(probe_x),  declines  when  the  probe  moves  away  from  the  Y-axis.  To 
explain  this  phenomena,  we  have  to  refer  back  to  the  lower  sphere  positional  uncertainty 
because  the  probe  positional  uncertainty  is  mainly  determined  by  the  uncertainty  of  the 
lower  sphere.  By  modifying  the  Equation  7-13,  the  lower  sphere  positional  uncertainty  in 
X-coordinate  can  be  written  as 


(I  \ 
\  r  ) 


(7,  + 


(I  \ 


( 


-1 1/2 


+  2 


(l  I  } 

l2l4 


\  r  J 


\r  J 
cov(A/2A/4)  +  2 


V 


2r2 


+  ■ 


v 


(I  \ 
L 

yr  J 


2rl 


cov(A/2Ar) 


'  J 


+  2 


\r  J 


-i.W-h  )  .  l 


2r2 


+  — 

2 


cov(A/4Ar) 


(7-20) 
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Probe  position  Uncertainty  -Uc(Probe_x) 


Probe  position  Uncertainty  -Uc(Probe_y) 


Probe  position  Uncertainty  -Uc(Probe_z) 


Figure  7- 1 0  The  combined  standard  uncertainty  for  probe  due 
to  the  thermal  error  and  strut  displacement  measurement  error 
( the  probe  at  the  level  of  Z=-25) 
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In  Equation  7-20,  the  terms  with  covariance  can  be  eliminated  since  there  is  no 
covariance  between  the  parameters  of  l2,  U,  and  r.  From  Table  7-5  we  know  that  cr^ 


=cf  and  <7r2= 1.79cr^ .  Then,  Equation  7-20  can  be  written  as 


(I  \ 


Gn  + 


(I  \ 

4 


<T,2+1.79 


V 


2r' 


W2 


V 


1/2 


(7-21) 


Without  the  last  term  in  Equation  7-21,  one  can  observe  that  the  combined 
standard  uncertainty  in  the  lower  sphere  in  X-coordinate  will  show  symmetry  to  the  plane 
of  X=34.1,  as  shown  in  Figure  7-11.  When  the  last  term  in  Equation  7-21  was 
considered,  we  observed  that  the  closer  the  probe  to  the  Y-axis,  the  larger  the  magnitude 
in  the  last  term.  This  is  the  reason  why  we  found  the  probe's  positional  uncertainty  in  the 
X-coordinate  declines  when  the  probe  moves  away  from  the  Y-axis. 


Figure  7-1 1  The  lower  sphere  positional  uncertainty  without 
the  last  term  in  Equation  7-21 
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Probe  Positional  Uncertainty  Due  to  Cylindrical  Column 

It  is  known  that  the  thermal  growth  in  the  cylindrical  column  will  be  directly 
added  onto  the  probe  measurement  error  in  the  Z-direction,  but  not  through  the  error 
propagation  procedure.  From  Table  6-5,  the  maximum  thermal  growth  in  the  column  is 
56.84  |iin  under  the  ambient  temperature  frequency  of  2  hour  per  cycle  and  the  variation 
amplitude  of  1 .4°F.  Therefore,  the  standard  uncertainty  of  the  probe  measurement  in  the 
Z-direction  due  to  the  thermal  growth  in  the  cylindrical  column  will  be  32.8  lfxin. 


In  the  previous  sections  we  investigated  the  probe  positional  uncertainty  due  to 
three  major  error  sources.  What  would  be  the  overall  combined  standard  uncertainty  for 
the  probe?  Here  we  will  discuss  two  extreme  cases: 

1 .  These  three  sources  of  uncertainty  are  independent  each  other. 

2.  These  three  sources  of  uncertainty  are  hundred  percent  correlated. 

For  the  first  case  the  overall  combined  standard  uncertainty  can  be  obtained  by 
taking  the  square  root  of  the  sum  of  the  squares  of  each  individual  standard  uncertainty 
(RSS)  as  shown  in  following  equation: 


For  the  second  case  the  combined  standard  uncertainty  for  the  probe  would  be  the 
sum  of  each  individual  standard  uncertainty: 


As  we  know,  the  overall  probe  combined  standard  uncertainty  must  stay  within 
the  two  boundaries.  Thus,  the  overall  combined  standard  uncertainty  for  the  probe 


7.5  The  Overall  Probe  Positional  Uncertainty 


(7-22) 


Uc  =  W,  +  U2  +  U 


(7-23) 
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position  is  determined  by  taking  the  average  of  RSS  and  the  individual  sum.  In  Table  7- 
7,  we  listed  the  worst  case  of  each  individual  uncertainty  sources.  From  Table  7-7  we 
observed  that  the  probe  Y-coordinate  has  the  smallest  positional  uncertainty  and  the  X 
coordinate  has  the  largest  uncertainty.  To  improve  the  probe's  positional  uncertainty,  one 
can  either  improve  the  hardware  so  that  calibration  uncertainty  could  be  reduced,  or  re- 
calibrate the  kinematic  parameter  when  the  thermal  growth  is  too  big. 


Table  7-7  The  uncertainty  budget  for  probe  position 


sources  of  uncertainty 

uc(Probe_x) 
(Win) 

uc(Probe_y) 
(Win) 

uc(Probe_z) 
(Win) 

1 

Calibration  (for  m=30) 

190.3 

69.5 

45 

2 

Thermal  growth  in  triangle 
frame  &  strut  linear 
displacement  measurement 

26.5 

24 

99.6 

3 

Thermal  growth  in  cylindrical 
column 

0 

0 

32.8 

Sum 

216.8 

93.5 

177.4 

RSS 

192.1 

73.5 

114.1 

Average(Sum+RSS) 

204.45 

83.5 

145.75 

It  has  to  be  noted  that  the  results  shown  in  Table  7-7  is  based  on  the  worst  case, 
operation  in  the  metrology  laboratory,  and  no  re-calibration  procedure.  To  obtain  the 
probe's  positional  uncertainty  for  different  ambient  environments,  one  can  follow  the 
procedures  we  have  performed  from  Sections  7-1  to  7-5. 


CHAPTER  8 
CONCLUSIONS  AND  FUTURE  WORK 

8.1  Conclusions 

A  self-calibration  methodology  for  a  hexapod  coordinate  measuring  machine  has 
been  developed.  It  was  found  that  the  self-calibration  technique  is  very  sensitive  to  the 
strut  linear  displacement  measurement  noise.  It  has  been  shown  that  the  calibration 
simulation  results  are  linearly  decreased  with  the  increase  of  the  strut  displacement 
measurement  uncertainty. 

To  minimize  the  effect  of  inevitable  measurement  noise  on  the  results  of 
estimation,  the  selection  of  measurement  configuration  becomes  very  important.  Thus, 
the  condition  number  was  used  to  be  the  index  to  quantify  the  goodness  of  pose  selection. 

A  procedure  to  determine  the  time  period  that  the  HCMM  needs  to  be  re- 
calibrated or  a  part  measurement  has  to  be  completed  has  been  demonstrated.  The 
combined  standard  uncertainty  analysis  has  been  performed  to  evaluate  the  HCMM's 
measurement  uncertainty. 

8.2  Future  Work 

Several  tasks  will  need  to  be  continually  performed  in  order  to  build  a  successful 
machine.  First,  the  methodology  for  verifying  the  probe  alignment  has  to  be  developed. 
This  will  avoid  the  error  source  come  from  the  probe  misalignment  from  the  center  rod. 
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Second,  the  methodology  for  verifying  the  HCMM  accuracy  needs  to  be  established. 
This  procedure  will  be  used  to  confirm  the  theoretical  analysis  results  as  shown  in 
Chapter  7.  Third,  improving  the  center  rod  design  is  necessary  so  that  the  thermal  growth 
in  the  center  rod  during  the  calibration  procedure  can  be  minimized.  Forth,  the  structural 
dynamic  analysis  has  to  be  performed  so  that  operation  of  the  machine  at  the  speed  of  its 
natural  frequency  can  be  avoided. 


APPENDIX  A 
OBSERVABILITY  INDICES 


Minimum  Singular  Values  and  Observability  Index 


For  a  closed  loop  kinematic  like  that  of  the  HCMM,  the  linearized  error  function 
can  be  written  as  e  -  Ja.  From  matrix  norm  basics  one  can  derive  the  following  ratio: 


where  ar  and  <7,  are  the  minimum  and  maximum  singular  of  matrix  J.  It  is  desired  that 
small  errors  in  a  will  make  a  largest  effect  on  the  residual  error  function  of  e.  Thus  we 
wish  to  make      /  |a|  as  big  as  possible.  This  implies  that  the  bigger  the  smallest 
singular  value  or  the  larger  the  sensitivity. 

From  a  geometrical  point  of  view,  if  we  assume  that  a  defines  a  hypersphere  with 
a  unit  radius,  then  e  is  a  hyperellipsoid  whose  semi-axes  are  the  singular  values  of  J. 
Thus  the  geometric  rationale  is  to  make  the  volume  of  the  hyperellipsoid  as  big  as 
possible.  Borm  and  Meng  [14]  prove  that  the  volume  is  proportional  to  the  value  as 
follows: 


obi  = 


where  obi  is  the  observability  index  and  m  is  the  number  of  poses. 
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Noise  amplification  Index 

This  number  is  the  combination  of  condition  number  and  minimum  singular  value 
which  is  given  as 

cr  cr2 

O  =  : —  =  — 

cond(J)  a, 

Nahvi  and  Hollerbach  [16]  claimed  that  this  index  is  more  sensitive  than  the  other 
three  indices. 


APPENDIX  B 
OPTIMAL  POSE  SETS 
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APPENDIX  C 
MATLAB  COMMAND  FILES 


The  Non-Linear  Least-Squares  Methods 


%%  1.  The  objective  is  to  calibrate  the  HCMM's  kinematic  parameters  using  a 
%%     linearized  least-squares  method. 

%%  2.  The  Jacobian  Matrix  [J]  was  calculated  by  using  the  analytical  deviation. 
%%  3.  The  input  information: 

%%        x:  the  initial  guess  for  kinematic  parameters;  size  (1x10) 
%%        m:  number  of  poses  used  in  calibration;  size  (1x1) 
%%        s:  the  strut  displacement  data;  size  (6xm) 
%%  4.  The  output  data: 

%%       a:  the  identified  kinematic  parameters;  size  (1x10) 
%%       J:  the  Jacobian  matrix;  size  (mx10) 


function  [a,J,jump]=hcmm_nls(x,m,s) 

count=1; 

n=length(x); 

%%%  r=a(7)  %%% 
%%%  b=a(8)  %%% 
%%%  h=a(9)  %%% 

%%%  d=a(10)  inital  guess  for  center  rod  length%%% 

repeat=1 ; 

stop=1; 

err_old=1e10; 

change_a5=0; 
change_a6=0; 
jump=0; 

adpt1=1;  adpt2=1; 

while  repeat  ==1  %%  preventing  a  negative  in  square  root  %% 
a=x; 

if  change_a5==1 

a(5)  =  a(5)+0.2*adpt1 ;  %increase  a5  to  prevent  negative  sqrt(-) 
elseif  change_a6==2 

a(6)  =  a(6)+0.2*adpt2; 

end 

change_a5=  0;  change_a6=0; 
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while  stop==1 ,  %%  stop  iteration  if  max(abs(dphi))  <  1e-12  %% 
if  jump  >  10  %%  automatically  adjust  the  initial  values  of  a5  &  a6  for  10  times 

repeat=0; 

stop=0; 

break 

end 

for  k=1:m; 

Xu(k)=  ( (a(1)+s(1,k))A2-(a(3)+s(3,k))A2+a(7)A2  )/(2*a(7)); 
Xl(k)=  ( (a(2)+s(2,k))A2-(a(4)+s(4,k))A2+a(7)A2  )/(2*a(7)); 

Yu(k)=  ( (a(3)+s(3,k))A2  -  (a(5)+s(5,k))A2  +2*Xu(k)*(a(7)-a(8))-(a(7)A2-a(8)A2-a(9)A2)) 
/(2*a(9)); 

Yl(k)=  ( (a(4)+s(4,k))A2  -  (a(6)+s(6,k))A2  +2*XI(k)*(a(7)-a(8))-(a(7)A2-a(8)A2-a(9)A2))/(2*a(9)); 
adj11=  (a(1)+s(1,k))A2  -  Xu(k)A2  -  Yu(k)A2  ; 
adj21=  (a(2)+s(2,k))A2  -  XI(k)A2  -  YI(k)A2; 
if  adj11  >0 

uuuuz=adj1 1 ; 
else 

test1=adj1 1 ;  %%  ->the  following  steps  are  used  to  automatically  adjust  the  initial  — 
change_a5=1 ;  %    values  to  prevent  the  negative  value  in  the  square  root  of  Z  I 
adpt1=adpt1+1; 
err_old=1e10; 
jump=jump+1; 
break 
end 

if  adj21  >  0 

Illlz=adj21; 
else 

test2=adj21 ; 

change_a6=2; 

adpt2=adpt2+1; 

err_old=1e10; 

jump=jump+1;%automatically  adjust  the  initial  value  if  negative  values  appear  in  adj1 1&2 

break  %       %  t 

end 

Zu(k)=  sqrt(  uuuuz ); 
Zl(k)=  -sqrt(llllz ); 

Lc(k)=  sqrt(  (Xu(k)-XI(k))A2  +  (Yu(k)-YI(k))A2  +  (Zu(k)-ZI(k))A2); 
F(k)=Lc(k)A2-a(10)A2; 
errorF(k)=  O.-F(k); 
end 

if  jump  >  10 
repeat=0; 
stop=0; 
break 

end 

if  change_a5  +change_a6  >  0 

break 
end 

err(count)=errorF*errorF';  %  Objective  value 
a_old=a; 

%%  Take  the  derivative  of  F 
for  k=1  :m 
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Xu_a1(k)=(a(1)+s(1,k))/a(7); 
Xu_a2(k)=0; 

Xu_a3(k)=-(a(3)+s(3,k))/a(7); 
Xu_a4(k)=0; 
Xu_a5(k)=0; 
Xu_a6(k)=0; 

Xu_r(k)=-0.5*((a(1  )+s(1  ,k))A2-(a(3)+s(3,k))A2)/a(7)A2+0.5; 

Xu_b(k)=0; 

Xu_h(k)=0; 

Xl_a1(k)=0; 

Xl_a2(k)=(a(2)+s(2,k))/a(7); 
Xl_a3(k)=0; 

Xl_a4(k)=-(a(4)+s(4,k))/a(7); 

Xl_a5(k)=0; 

Xl_a6(k)=0; 

XI_r(k)=-0.5*((a(2)+s(2,k))A2-(a(4)+s(4,k))A2)/a(7)A2+0.5; 

Xl_b(k)=0; 

Xl_h(k)=0; 

Yu_a1  (k)=(a(7)-a(8))/a(9)*(Xu_a1  (k)); 
Yu_a2(k)=0; 

Yu_a3(k)=-a(8)/a(9)*Xu_a3(k); 
Yu_a4(k)=0; 

Yu_a5(k)=-(a(5)+s(5,  k))/a(9) ; 
Yu_a6(k)=0; 

Yu_r(k)=Xu(k)/a(9)-  a(7)/a(9)  +(a(7)-a(8))/a(9)*(Xu_r(k)); 
Yu_b(k)=-Xu(k)/a(9)+a(8)/a(9); 

Yu_h(k)=-0.5*(  (a(3)+s(3,k))A2  -  (a(5)+s(5,k))A2  +2*Xu(k)*(a(7)-a(8)) 
-(a(7)A2-a(8)A2)  )/a(9)A2+0.5; 

Yl_a1(k)=0; 

YI_a2(k)=(a(7)-a(8))/a(9)*(XI_a2(k)); 
Yl_a3(k)=0; 

YI_a4(k)=-a(8)/a(9)*XI_a4(k); 
Yl_a5(k)=0; 

Yl_a6(k)=-(a(6)+s(6,k))/a(9); 

YI_r(k)=XI(k)/a(9)  -  a(7)/a(9)  +(a(7)-a(8))/a(9)*(XI_r(k)); 
YI_b(k)=-XI(k)/a(9)+a(8)/a(9); 

YI_h(k)=-0.5*(  (a(4)+s(4,k))A2  -  (a(6)+s(6,k))A2  +2*XI(k)*(a(7)-a(8)) 
-(a(7)A2-a(8)A2)  )/a(9)A2+0.5; 

asxy(k)=1/sqrt((a(1)+s(1,k))A2-Xu(k)A2-Yu(k)A2); 

Zu_a1  (k)=((a(1  )+s(1  ,k))-Xu(k)*Xu_a1  (k)-Yu(k)*Yu_a1  (k))*asxy(k); 

Zu_a2(k)=0; 

Zu_a3(k)=(-Xu(k)*Xu_a3(k)-Yu(k)*Yu_a3(k))*asxy(k); 
Zu_a4(k)=0; 

Zu_a5(k)=(-Yu(k)*Yu_a5(k))*asxy(k); 
Zu_a6(k)=0; 

Zu_r(k)=(-Xu(k)*Xu_r(k)-Yu(k)*Yu_r(k))*asxy(k); 

Zu_b(k)=(-Yu(k)*Yu_b(k))*asxy(k); 

Zu_h(k)=(-Yu(k)*Yu_h(k))*asxy(k); 

asxYI(k)=-1/sqrt((a(2)+s(2,k))A2-XI(k)A2-YI(k)A2); 
Zl_a1(k)=0; 
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ZI_a2(k)=((a(2)+s(2,k))-XI(k)*XI_a2(k)-YI(k)*YI_a2(k))*asxYI(k); 
Zl_a3(k)=0; 

ZI_a4(k)=(-XI(k)*XI_a4(k)-YI(k)*YI_a4(k))*asxYI(k); 
Zl_a5(k)=0; 

ZI_a6(k)=(-YI(k)*YI_a6(k))*asxYI(k); 
ZI_r(k)=(-XI(k)*XI_r(k)-YI(k)*YI_r(k))*asxYI(k); 
ZI_b(k)=(-YI(k)*YI_b(k))*asxYI(k); 
ZI_h(k)=(-YI(k)*YI_h(k))*asxYI(k); 

xyz(k)=1.0;  %  this  term  is  for  objective  function  F=LceA2-a(1 0)A2 
Xul(k)=Xu(k)-XI(k);  Yul(k)=Yu(k)-YI(k);  Zul(k)=Zu(k)-ZI(k); 

Lce_a1  (k)=(Xul(k)*(Xu_a1  (k)-XI_a1  (k))+Yul(k)*(Yu_a1  (k) 

-Yl_a1  (k))+Zul(k)*(Zu_a1  (k)-ZI_a1  (k)))*xyz(k); 
Lce_a2(k)=(Xul(k)*(Xu_a2(k)-XI_a2(k))+Yul(k)*(Yu_a2(k) 

-YI_a2(k))+Zul(k)*(Zu_a2(k)-ZI_a2(k)))*xyz(k); 
Lce_a3(k)=(Xul(k)*(Xu_a3(k)-XI_a3(k))+Yul(k)*(Yu_a3(k) 

-YI_a3(k))+Zul(k)*(Zu_a3(k)-ZI_a3(k)))*xyz(k); 
Lce_a4(k)=(Xul(k)*(Xu_a4(k)-XI_a4(k))+Yul(k)*(Yu_a4(k) 

-YI_a4(k))+Zul(k)*(Zu_a4(k)-ZLa4(k)))*xyz(k); 
Lce_a5(k)=(Xul(k)*(Xu_a5(k)-XI_a5(k))+Yul(k)*(Yu_a5(k) 

-YI_a5(k))+Zul(k)*(Zu_a5(k)-ZI_a5(k)))*xyz(k); 
Lce_a6(k)=(Xul(k)*(Xu_a6(k)-XI_a6(k))+Yul(k)*(Yu_a6(k) 

-YI_a6(k))+Zul(k)*(Zu_a6(k)-ZI_a6(k)))*xyz(k); 
Lce_r(k)-(Xul(k)*(Xu_r(k)-XI_r(k))+Yul(k)*(Yu_r(k)-YI_r(k)) 

+Zul(k)*(Zu_r(k)-ZI_r(k)))*xyz(k); 
Lce_b(k)=(Xul(k)*(Xu_b(k)-XLb(k))+Yul(k)*(Yu_b(k)-YI_b(k)) 

+Zul(k)*(Zu_b(k)-ZI_b(k)))*xyz(k); 
Lce_h(k)=(Xul(k)*(Xu_h(k)-XI_h(k))+Yul(k)*(Yu_h(k)-YI_h(k)) 

+Zul(k)*(Zu_h(k)-ZI_h(k)))*xyz(k); 

f_a1(k)=2*Lce_a1(k); 

f_a2(k)=2*Lce_a2(k); 

f_a3(k)=2*Lce_a3(k); 

f_a4(k)=2*Lce_a4(k); 

f_a5(k)=2*Lce_a5(k); 

f_a6(k)=2*Lce_a6(k); 

f_r(k)=2*Lce_r(k); 

f_b(k)=2*Lce_b(k); 

f_h(k)=2*Lce_h(k); 

f_Lct(k)=-2*a(10);  %this  term  is  for  bjective  function  F=LceA2-a(1 0)A2 


J=[La1\  f_a2',  f_a3',  f_a4',  f_a5',  f_a6',  f_r\  f_b',  f_h',  f_Lct'];  %  the  Jacobian  matrix 
%%%  end  deriviation  %%%% 
k1=cond(J);  %  condition  number 

if  (k1)  >  1e5  %  if  condition  number  >  1e5  ,  then  terminate  program 

disp('  condition  number  too  big.  (ill  condition)'); 

disp('  press  CTR-C  to  terminate  program'); 
end 

%%%%  find  increment  %%%% 
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dphi=inv(J'\J)\TerrorF'; 


a=a_old  +  dphi';  %%%  update  initial  parameters  of  a  %%% 
max(dphi); 

if  max(abs(dphi))  <  1e-12  %%%  terminate  iteration 
q=9999999999999 
stop=0; 
repeat=0; 

elseif  err_old  <  err(count)  %%  terminate  iteration  if  the  previous  objective  value  is  less  than 

q=88888888888888     %%  the  current  one. 

stop=0; 

repeat=0; 
end 

err_old=err(count); 
x=a; 

count=count+1 ;  %%%  count  the  iteration  %%% 
end 
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The  Probe  Positional  Uncertainty  Evaluation  File 


%  File  name:  "hcm_unct.m" 

%  This  program  is  to  evaluate  the  HCMM's  probe  positional  uncertainty 

%  Input  files: 

%      verif_p.dat:  The  probe  positions  where  we  want  to  evaluate  the  uncertainty.  Matrix  size 
%  (m  x  6) 

%      verif_dl.dat:  The  strut  displacement  data  at  each  probe  position.  Matrix  size:(m  x6) 
%      covar_mtx.dat:  The  covariance  between  each  kinematic  parameters.  Matrix  size:(9  x  9) 

clear 

load  verif_p.dat;  %  load  the  positional  coordinates  for  evaluation;(Xu,  Yu,  Zu,  XI,  Yl,  Zl);size  (mx6) 
p1  p2=verif_p';    %  transpose  into  (6xm) 

load  verif_dl.dat;  %  load  the  strut  displacement  data;  (dl_1,  dl_2,  dl_3,  dl_4,  dl_5,  dL6);sixe(mx6) 
s=verif_dl';         %  transpose  the  matrix  into  (6xm) 

load  covar_mtx.dat;  %  load  the  covariance  between  each  parameters;  see  Table7-3  &7-5 
cov_a=covar_mtx; 

X(1 :6)=43.10642755243507*ones(1 ,6); 
X(7)=68.233;  X(8)=X(7)/2;  X(9)=X(8)*sqrt(3); 


M=size(p1p2,2); 
a=X; 

Xu=p1  p2(1  ,:);Yu=p1  P2(2,:);Zu=p1  p2(3,:);XI=p1  P2(4,:);YI=p1  p2(5,:);ZI=p1  p2(6. :); 
VarXu_old=0; 


for  k=1  :M 

Xu_a1(k)=(a(1)+s(1,k))/a(7); 
Xu_a2(k)=0; 

Xu_a3(k)=-(a(3)+s(3,k))/a(7); 
Xu_a4(k)=0; 
Xu_a5(k)=0; 
Xu_a6(k)=0; 

Xu_r(k)=-0.5*((a(1  )+s(1  ,k))A2-(a(3)+s(3,k))A2)/a(7)A2+0.5; 

Xu_b(k)=0; 

Xu_h(k)=0; 

Xl_a1(k)=0; 

XLa2(k)=(a(2)+s(2,k))/a(7); 
Xl_a3(k)=0; 

Xl_a4(k)=-(a(4)+s(4,k))/a(7); 

Xl_a5(k)=0; 

Xl_a6(k)=0; 

XI_r(k)=-0.5*((a(2)+s(2,k))A2-(a(4)+s(4Ik))A2)/a(7)A2+0.5; 

Xl_b(k)=0; 

Xl_h(k)=0; 

Yu_a1  (k)=(a(7)-a(8))/a(9)*(Xu_a1  (k)); 
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Yu_a2(k)=0; 

Yu_a3(k)=-a(8)/a(9)*Xu_a3(k); 
Yu_a4(k)=0; 

Yu_a5(k)=-(a(5)+s(5,k))/a(9); 
Yu_a6(k)=0; 

Yu_r(k)=Xu(k)/a(9)-  a(7)/a(9)  +(a(7)-a(8))/a(9)*(Xu_r(k)); 
Yu_b(k)=-Xu  (k)/a(9)+a(8)/a(9) ; 

Yu_h(k)=-0.5*(  (a(3)+s(3,k))A2  -  (a(5)+s(5,k))A2  +2*Xu(k)*(a(7)-a(8)Ha(7)A2-a(8)A2) ) 
/a(9)A2+0.5; 

Yl_a1(k)=0; 

YLa2(k)=(a(7)-a(8))/a(9)*(XI_a2(k)); 
Yl_a3(k)=0; 

YI_a4(k)=-a(8)/a(9)*XI_a4(k); 
Yl_a5(k)=0; 

Yl_a6(k)=-(a(6)+s(6,k))/a(9); 

YI_r(k)=XI(k)/a(9)  -  a(7)/a(9)  +(a(7)-a(8))/a(9)*(XI_r(k)); 
YI_b(k)=-XI(k)/a(9)+a(8)/a(9); 

YI_h(k)=-0.5*(  (a(4)+s(4,k))A2  -  (a(6)+s(6,k))A2  +2*XI(k)*(a(7)-a(8))-(a(7)A2-a(8)A2) ) 
/a(9)A2+0.5; 

asxy(k)=1/sqrt((a(1)+s(1,k))A2-Xu(k)A2-Yu(k)A2); 

Zu_a1  (k)=((a(1  )+s(1  ,k))-Xu(k)*Xu_a1  (k)-Yu(k)*Yu_a1  (k))*asxy(k); 

Zu_a2(k)=0; 

Zu_a3(k)=(-Xu(k)*Xu_a3(k)-Yu(k)*Yu_a3(k))*asxy(k); 
Zu_a4(k)=0; 

Zu_a5(k)=(-Yu(k)*Yu_a5(k))*asxy(k); 
Zu_a6(k)=0; 

Zu_r(k)=(-Xu(k)*Xu_r(k)-Yu(k)*Yu_r(k))*asxy(k); 

Zu_b(k)=(-Yu(k)*Yu_b(k))*asxy(k); 

Zu_h(k)=(-Yu(k)*Yu_h(k))*asxy(k); 

asxYI(k)=-1/sqrt((a(2)+s(2,k))A2-XI(k)A2-YI(k)A2); 
Zl_a1(k)=0; 

ZI_a2(k)=((a(2)+s(2,k))-XI(k)*XI_a2(k)-YI(k)*YLa2(k))*asxYI(k); 
Zl_a3(k)=0; 

ZI_a4(k)=(-XI(k)*XI_a4(k)-YI(k)*YLa4(k))*asxYI(k); 
Zl_a5(k)=0; 

ZI_a6(k)=(-YI(k)*YI_a6(k))*asxYI(k); 
ZI_r(k)=(-XI(k)*XI_r(k)-YI(k)*YI_r(k))*asxYI(k); 
ZI_b(k)=(-YI(k)*YI_b(k))*asxYI(k); 
ZI_h(k)=(-YI(k)*YI_h(k))*asxYI(k); 

xyz(k)=1/sqrt((Xu(k)-XI(k))A2+(Yu(k)-YI(k))A2+(Zu(k)-ZI(k))A2); 
Xul(k)=Xu(k)-XI(k);  Yul(k)=Yu(k)-YI(k);  Zul(k)=Zu(k)-ZI(k); 

Lc_a1  (k)=(Xul(k)*(Xu_a1  (k)-XI_a1  (k))+Yul(k)*(Yu_a1  (k)-YI_a1  (k))+Zul(k)*(Zu_a1  (k) 
-Zl_a1(k)))*xyz(k); 

Lc_a2(k)=(Xul(k)*(Xu_a2(k)-XI_a2(k))+Yul(k)*(Yu_a2(k)-YI_a2(k))+Zul(k)*(Zu_a2(k) 
-Zl_a2(k)))*xyz(k); 

Lc_a3(k)=(Xul(k)*(Xu_a3(k)-XI_a3(k))+Yul(k)*(Yu_a3(k)-YI_a3(k))+Zul(k)*(Zu_a3(k) 
-Zl_a3(k)))*xyz(k); 

Lc_a4(k)=(Xul(k)*(Xu_a4(k)-XI_a4(k))+Yul(k)*(Yu_a4(k)-YI_a4(k))+Zul(k)*(Zu_a4(k) 
-Zl_a4(k)))*xyz(k); 

Lc_a5(k)=(Xul(k)*(Xu_a5(k)-XI_a5(k))+Yul(k)*(Yu_a5(k)-YI_a5(k))+Zul(k)*(Zu_a5(k) 
-Zl_a5(k)))*xyz(k); 

Lc_a6(k)=(Xul(k)*(Xu_a6(k)-XI_a6(k))+Yul(k)*(Yu_a6(k)-YI_a6(k))+Zul(k)*(Zu_a6(k) 
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-Zl_a6(k)))*xyz(k); 

Lc_r(k)=(Xul(k)*(Xu_r(k)-XI_r(k))+Yul(k)*(Yu_r(k)-YI_r(k))+Zul(k)*(Zu_r(k) 
-Zl_r(k)))*xyz(k); 

Lc_b(k)=(Xul(k)*(Xu_b(k)-XI_b(k))+Yul(k)*(Yu_b(k)-YI_b(k))+Zul(k)*(Zu_b(k) 
-Zl_b(k)))*xyz(k); 

Lc_h(k)=(Xul(k)*(Xu_h(k)-XI_h(k))+Yul(k)*(Yu_h(k)-YI_h(k))+Zul(k)*(Zu_h(k) 
-Zl_h(k)))*xyz(k); 

end 

fork=1:M 
VarXu_old=0; 
VarYu_old=0; 
VarZu_old=0; 
VarXI_old=0; 
VarYI_old=0; 
VarZI_old=0; 
VarLc_old=0; 

U=[Xu_a1(k)  Xu_a2(k)  Xu_a3(k)  Xu_a4(k)  Xu_a5(k)  Xu_a6(k)  Xu_r(k)  Xu_b(k)  Xu_h(k) 
Yu_a1(k)  Yu_a2(k)  Yu_a3(k)  Yu_a4(k)  Yu_a5(k)  Yu_a6(k)  Yu_r(k)  Yu_b(k)  Yu_h(k) 
Zu_a1(k)  Zu_a2(k)  Zu_a3(k)  Zu_a4(k)  Zu_a5(k)  Zu_a6(k)  Zu_r(k)  Zu_b(k)  Zu_h(k) 
Xl_a1(k)  Xl_a2(k)  Xl_a3(k)  Xl_a4(k)  Xl_a5(k)  Xl_a6(k)  Xl_r(k)  Xl_b(k)  Xl_h(k) 
Yl_a1(k)  Yl_a2(k)  Yl_a3(k)  Yl_a4(k)  Yl_a5(k)  Yl_a6(k)  Yl_r(k)  Yl_b(k)  Yl_h(k) 
Zl_a1(k)  Zl_a2(k)  Zl_a3(k)  Zl_a4(k)  Zl_a5(k)  Zl_a6(k)  Zl_r(k)  Zl_b(k)  Zl_h(k) 
Lc_a1  (k)  Lc_a2(k)  Lc_a3(k)  Lc_a4(k)  Lc_a5(k)  Lc_a6(k)  Lc_r(k)  Lc_b(k)  Lc_h(k)]; 

for  i=1:size(U,2)  %%  see  Equation  7-9  &  7-10  for  calculating  the  variances 
forj=1:size(U,2) 

VarXu(k)=U(1  ,i)*U(1  J)*cov_a(i,j)+VarXu_old; 

VarYu(k)=U(2,i)*U(2,j)*cov_a(i,j)+VarYu_old; 

VarZu(k)=U(3,i)*U(3,j)*cov_a(i,j)+VarZu_old; 

VarXu_old=VarXu(k); 

VarYu_old=VarYu(k); 

VarZu_old=VarZu(k); 

VarXI(k)=U(4,i)*U(4,j)*cov_a(i,j)+VarXI_old; 

VarYI(k)=U(5,i)*U(5,j)*cov_a(i,j)+VarYI_old; 

VarZI(k)=U(6,i)*U(6Ij)*cov_a(i,j)+VarZI_old; 

VarXI_old=VarXI(k); 

VarYI_old=VarYI(k); 

VarZI_old=VarZI(k); 

VarLc(k)=U(7,i)*U(7,j)*cov_a(i,j)+VarLc_old; 

VarLc_old=  VarLc(k); 
end 
end 

VarpobsumX_old=0; 
VarpobsumY_old=0; 
VarpobsumZ_old=0; 
for  i=1:size(U,2) 
forj=1:size(U,2) 

VarpobsumX=U(1,i)*U(4,j)*cov_a(i,j)+  VarpobsumX_old;  %last  term  of  Eqn.7-19 
VarpobsumY=U(2,i)*U(5,j)*cov_a(i,j)+  VarpobsumY_old; 
VarpobsumZ=U(3,i)*U(6,j)*cov_a(i,j)+  VarpobsumZ_old; 
VarpobsumX_old=VarpobsumX; 
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VarpobsumY_old=VarpobsumY; 
VarpobsumZ_old=VarpobsumZ; 
end 
end 

Var_probe_X(k)=(1+5/35)A2*VarXI(k)+(5/35)A2*VarXu(k)-2*(1+5/35)*(5/35)*VarpobsumX; 
Var_probe_Y(k)=(1+5/35)A2*VarYI(k)+(5/35)A2*VarYu(k)-2*(1+5/35)*(5/35)*VarpobsumY; 
Var_probe_Z(k)=(1+5/35)A2*VarZI(k)+(5/35)A2*VarZu(k)-2*(1+5/35)*(5/35)*VarpobsumZ; 


StdXu=sqrt(VarXu); 

StdYu=sqrt(VarYu); 

StdZu=sqrt(VarZu); 

StdXI=sqrt(VarXI); 

StdYI=sqrt(VarYI); 

StdZI=sqrt(VarZI); 

StdLc=sqrt(VarLc); 

Std  P  robe_x=sq  rt  (Va  r_p  robe_X) ; 
StdProbe_y=sqrt(Var_probe_Y); 
StdProbe_z=sqrt(Var_probe_Z); 

for  i=1:M 

dPu(i)=sqrt(StdXu(i)A2+StdYu(i)A2+StdZu(i)A2); 
dPI(i)=sqrt(StdXI(i)A2+StdYI(i)A2+StdZI(i)A2); 
dProbe(i)=sqrt(StdProbe_x(i)A2+StdProbe_y(i)A2+StdProbe_z(i)A2); 
end 

%  The  following  codes  are  used  to  plot  the  positional  uncertainties  graphics 
k=1; 

for  i=1:size(p1p2,1) 
forj=1:size(p1p2,2) 

dd1(j,i)=StdXu(k); 

dd2(j,i)=StdYu(k); 

dd3(j,i)=StdZu(k); 

dd4(j,i)=StdLc(k); 

dd5(j,i)=StdXI(k); 

dd6(j,i)=StdYI(k); 

dd7(j,i)=StdZI(k); 

dd8(j,i)=StdProbe_x(k); 

dd9(j,i)=StdProbe_y(k); 

dd10(j,i)=StdProbe_z(k); 

k=k+1; 
end 
end 

for  i=1:21  ^change  this  value  if  the  grid  in  Fig7-9  &  7-10  is  not  (21x21) 

XXu(i)=Xu(20*i+1); 
end 

for  i=1 :21  >  change  this  value  if  the  grid  in  Fig7-9  &  7-1 0  is  not  (21  x21 ) 

YYu(i)=Yu(i); 
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figure(1) 

surfc(XXu,YYu,dd1) 

xlabel('X') 

ylabel('Y') 

zlabel('Uc(Xu)  inch') 

title('Upper  Sphere  Position  Uncertainty  --Uc(Xu)') 
figure(2) 

surfc(XXu,YYu,dd2) 

xlabel('X') 

ylabel('Y') 

zlabel('Uc(Yu)  inch') 

title('Upper  Sphere  Position  Uncertainty  -Uc(Yu)') 
figure(3) 

surfc(XXu,YYu,dd3) 

xlabel('X') 

ylabel('Y') 

zlabel('Uc(Zu)  inch') 

title('Upper  Sphere  Position  Uncertainty  --Uc(Zu)') 
figure(4) 

surfc(XXu,YYu,dd4) 

xlabel('X') 

ylabel('Y') 

zlabel('Uc(Lc)  inch') 

%text(35,20,2.4e-5,,m=100  noise=+/-10e-6  inch') 
title('Center  Rod  Length  Uncertainty  -Uc(Lc)') 

figure(5) 

surfc(XXu,YYu,dd5) 

xlabel('X') 

ylabel('Y') 

zlabel('Uc(XI)  inch') 

title('Lower  Sphere  Position  Uncertainty  -Uc(XI)') 
figure(6) 

surfc(XXu,YYu,dd6) 

xlabel('X') 

ylabel('Y') 

zlabel('Uc(YI)  inch') 

title('Lower  Sphere  Position  Uncertainty  --Uc(YI)') 
figure(7) 

surfc(XXu,YYu,dd7) 

xlabeK'X') 

ylabel('Y') 

zlabel('Uc(ZI)  inch') 

title('Lower  Sphere  Position  Uncertainty  -Uc(ZI)') 
figure(8) 

surfc(XXu,YYu,dd8) 

xlabel('X') 

ylabel('Y') 

zlabel('Uc(Probe_x)  inch') 

title('Probe  Position  Uncertainty  --Uc(Probe_x)') 

figure(9) 

surfc(XXu,YYu,dd9) 

xlabel('X') 

ylabel('Y') 

zlabel('Uc(Probe_y)  inch1) 
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title('Probe  Position  Uncertainty  --Uc(Probe_y)') 
figure(10) 

surfc(XXu,YYu,dd10) 

xlabel('X') 

ylabel('Y') 

zlabel('Uc(Probe_z)  inch') 

title('Probe  position  Uncertainty  --Uc(Probe_z)') 
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